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The  objectives  of  this  woik  were  to  develop,  implement  and  validate  a  reaction  zone  model  for 
vorticity  based  turbulent  combustion  simulation  at  high  Reynolds  and  Damkohler  numbers.  Direct 
Simulation  results  using  the  transport  element  method  were  used  to  examine  the  structure  of  the 
reaction  zone  and  to  develop  a  reasonable  set  of  approximations  that  could  be  used  to  simplify  the 
governing  equations.  The  resulting  model,  adopting  a  singular  expansion  philosophy  of  the  flow 
equations;  the  elemental  flame  model  consists  of  (1)  a  conserved  scalar  approximation  of  the  outer 
non-reacting  flow  to  determine  the  location  of  the  reaction  surface;  (2)  an  unsteady,  uniformly 
strained  flame  structure  model  for  the  inner  reacting  flow  imbedded  within  the  reaction  surface,  to 
compute  the  local  burning  rate  and  flame  structure  profiles;  and,  (3)  a  set  of  kinematically  based 
approximations  used  to  monitor  the  generation,  interaction  and  elimination  of  flame  surface  area  as  it 
spins  around  and  reaches  the  tip  of  the  spiral  within  large  vortical  structures.  Comparisons  between 
these  “large  structure  simulations”  and  direct  numerical  simulation  showed  that  the  model  could 
accurately  capture  the  physics  of  reacting  flows  and  predict  flame  surface  evolution  and  rate  of 
burning.  Future  work  should  be  concerned  with  (a)  extending  this  model  by  incorporating  another 
approximation  at  areas  of  low  strains,  i.e.,  inside  the  large  stracture  where  the  reaction  zones  resemble 
those  of  stratified  reactors,  and  (b)  extending  the  application  of  the  developed  model  to  three- 

dimensional  flows. _ _  _  n  ' 
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ABSTRACT 

The  objectives  of  this  work  were  to  develop,  implement  and  validate  a  reaction  zone  model 
for  vorticity  based  turbulent  combustion  simulation  at  high  Reynolds  and  Damkohler 
numbers.  Direct  Simulation  results  using  the  transport  element  method  were  used  to 
examine  the  structure  of  the  reaction  zone  and  to  develop  a  reasonable  set  of 
approximations  that  could  be  used  to  simplify  the  governing  equations.  The  resulting 
model,  adopting  a  singular  expansion  philosophy  of  the  flow  equations;  the  elemental  flame 
model  consists  of  (1)  a  conserved  scalar  approximation  of  the  outer  non-reacting  flow  to 
determine  the  location  of  the  reaction  surface;  (2)  an  unsteady,  uniformly  strained  flame 
structure  model  for  the  inner  reacting  flow  imb^ded  within  Ae  reaction  surface,  to 
compute  the  local  burning  rate  and  flame  stracture  profiles;  and,  (3)  a  set  of  kinematically 
based  approximations  us^  to  monitor  the  generation,  interaction  and  elimination  of  flame 
surface  area  as  it  spins  around  and  reaches  the  tip  of  the  spiral  within  large  vortical 
structures.  Comparisons  between  these  “large  stracture  simulations”  and  direct  numerical 
simulation  showed  that  the  model  could  accurately  capture  the  physics  of  reacting  flows 
and  predict  flame  surface  evolution  and  rate  of  burning.  Future  work  should  be  concerned 
with  (a)  extending  this  model  by  incorporating  another  approximation  at  areas  of  low 
strains,  i.e.,  inside  the  large  stracture  where  the  reaction  zones  resemble  those  of  stratified 
reactors,  and  (b)  extending  the  application  of  the  developed  model  to  three-dimensional 
flows. 


«l  « 

I  TECHNICAL  REPORT 

1.1.  Obiectives 

The  objectives  of  this  project  were  to  develop,  implement  and  validate 
comprehensive  reaction  zone  models,  incorporating  chemical  kinetics  and  transport 
processes  at  a  level  of  detail  sufficient  for  the  prediction  of  nonequilibrium  combustion 
phenomena,  which  are  fully  compatible  with  Lagrangian  vortex  simulation  approaches, 
i.e.,  kinematically  based  flow  computations.  The  developed  models  would  be 
comprehensive  enough  to  enable  the  computation  of  in-flame  formation  of  minor  species  as 
well  as  post  flame  destruction/formation  of  these  species.  The  reaction  zone  models  would 
be  developed  for  the  high  Reynolds  number,  high  Damkohler  number  conditions,  i.e.,  in 
the  so-called  flamelet  regime,  or  what  we  call  the  elemental  flame  regime,  and  would  be 
compatible  with  large  eddy  simulation,  or  preferable  large  structure  simulations  of  turbulent 
reacting  flow.  The  models  would  be  guided  by  the  results  of  direct  numerical  simulation 
results  obtained  via  detailed  solutions  of  the  governing  equations  utilizing  the  transport 
element  method,  and  would  be  designed  to  converge  to  the  direct  simulation  results  as  the 
numerical  or  modeling  resolution  was  refined. 

1.2.  Overview 

We  have  accomplished  the  above  objectives  by  achieving  substantial  progress  in 
three  related  projects  summarized  next. 

(1)  Direct  Numerical  Simulation  of  Turbulent  combustion  at  High  Reynolds  Number  and 
Damkohler  Number  Using  the  T ransport  Element  Method. 

In  this  project,  we  showed  that  our  direct  simulation  method  can  be  extended  to 
accommodate  finite  rate  reaction  mechanisms  by  using  the  concept  of  product  elements  and 
by  refining  the  diffusion-mixing-reaction  algorithm  (Appendices  I  and  FV  and  Pubs.  3  and 
4).  Next  we  used  the  results  of  the  direct  simulations  to;  (a)  show  that  at  high  Damkohler 
number,  results  of  the  infinite  reaction  rate  model  using  Shvab-Zeldovich  (conserved 
scalar)  formulation  agree  well  with  those  obtained  using  a  finite  reaction  rate  model  in  terms 
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of  the  flow  dynamics,  the  mixing  rates  and  the  product  formation  rate,  (b)  gain  valuable 
insight  into  the  physics  of  turbulent  reacting  flows  to  guide  the  formxilation  of  the  reaction 
zone  models,  and  (c)  validate  the  results  obtained  using  the  reaction  zone  models  which  are; 
described  next. 

(2)  Reaction  Zone  Models  in  Lagrangian  Simulations  Using  the  Concept  of  Elemental 
Flames  at  Infinite  and  Finite  Reaction  Rates. 

In  this  project,  and  recognizing  that  direct  simulation  using  finite  rate  multispecies 
mechanisms  is  extremely  expensive  when  detailed  transport  processes  and  chemical 
kinetics  are  included  in  the  calculations,  we  embarked  on  an  effort  to  use  what  we  learned 
from  the  detailed  simulation  results  to  construct  physically  accurate  reaction  zone  models 
.  which  are  compatible  with  the  vortex  simulation  approach  of  solving  the  Navier-Stokes 
equation  and  are  applicable  at  high  Reynolds  and  Damkohler  numbers.  To  this  end: 

(a)  At  the  fast  chemistry  model,  i.e.  to  model  mixing  without  reaction,  we  formulated, 
implemented  and  validated  a  version  of  our  transport  element  method  which  takes 
advantage  of  the  observation  that  as  soon  as  the  mixing  layer  experiences  folding  to  the 
point  where  the  distance  between  neighboring  elements  is  smaller  than  the  diffusion 
thickness,  the  layers  cannot  be  distinguished  from  one  another  and  they  should  be 
merged.  Moreover,  as  soon  as  the  elements  reach  the  tip  of  the  rollup  spiral  inside  the 
large  eddies,  they  join  what  can  be  considered  a  well  mixed  zone,  or  a  macro  element  (it 
is  actually  a  radially  stratified  zone  but  with  small  gradients).  While  the  strain  and 
elongation  of  the  transport  elements  increase  their  number  as  they  move  around  the 
spiral  centers,  their  number  is  limited  by  elimination  at  the  spiral  tip. 

(b)  At  finite  rate  reactions,  instead  of  using  a  transport  element  which  simulates  the  mixing 
process  locally,  we  developed  a  new  concept:  an  elemental  flame.  This  is  a  reacting 
element  (molecule)  which  simulate  combustion  locally.  Elemental  flames  are 
distributed  along  the  spiraling  interface  to  capture  the  reaction  zone  normal  to  the 
interface.  Contrary  to  a  transport  element,  whose  field  is  obtained  analytically  as  the 
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Green  function  of  the  differential  equation  governing  convective-diffusive  transport,  the 
field  of  the  elemental  flame  is  obtained  numerically  by  integrating  a  reduced  set  of  the 
original  convection-diffusion-reaction  governing  equations  obtained  by  enforcing  a 
constant  value  of  the  strain  field  throughout  the  flame  structure  (while  these  equations 
represent  the  fitrst  term  in  a  series  expansion  in  the  strain,  the  effective  strain  used  in  the 
solution,  to  improve  accuracy,  is  selected  to  be  that  at  reaction-diffusion  zone,  as 
described  in  detail  in  Appendices  n  and  HI).  For  elemental  flames,  as  for  transport 
elements,  they  are  bom  at  the  zones  of  maximum  strain  (outside  the  large  eddies)  while 
they  die  where  the  strain  reaches  zero  (at  the  centers  of  the  eddies). 

(c)  Consistent  with  the  concept  of  elemental  flames  which  produce  macro  combustion 
elements  at  the  centers  of  the  eddies,  we  modified  the  vortex  simulation  scheme  to  take 
advantage  of  the  natural  tendency  of  vorticity  to  accumulate  inside  the  large  eddies  by 
replacing  a  collection  of  small  vortices  by  a  growing  large  element  This  “large 
structure  simulation”  scheme  has  been  shown  to  agree  very  well  with  the  results  of  the 
detailed  simulation  but  at  a  fi'action  of  the  cost 
(3)  Unsteady  Uniformly  Strained  Flame  Structure  Model  with  Detailed  Chemistry  and 
Transport  for  Computing  the  Elemental  Flame  Structure. 

As  mentioned  before  the  elemental  flame  model  of  the  reaction  zone  stmcture 
assumes  that  locally,  a  thin  flat  flame  exists.  Thus,  one  must  start  the  computation  by 
dividing  the  material  interface,  of  flame  into  a  large  number  of  flat  segments,  and  solve  a 
large  number  of  flame  structure  problems  to  reconstmct  the  reacting  flow  field.  To  do  this 
economically,  we  simplified  the  set  of  equations  governing  the  unsteady  strained  flame 
stmcture  by  assuming  that  the  strain  is  uniform  throughout  the  flame  stmcture  and  thus 
transforming  the  governing  equations  into  a  simpler  set  of  reaction-  diffusion  equations. 

By  properly  selecting  the  equivalent  uniform  strain  corresponding  to  the  given  actual  strain 
using  a  formula  which  we  have  derived  for  both  diffusion  and  premixed  flames,  one  can 
use  this  model,  which  is  more  robust  numerically  and  faster  computationally,  to  recover  the 
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flame  structure  accurately.  The  model  was  also  tested  and  validated  using  reduced  chemical 
kinetics  schemes. 

In  the  rest  of  the  report,  more  detail  summarizing  the  effort  leading  to  the 
formulation  of  the  elemental  flame  reaction  zone  models  are  presented.  The  Appendices 
describe  more  extended  versions  of  different  parts  of  the  projects. 

1.3.  Direct  Simulations  Approach  and  Results 

Direct  simulation  of  the  vorticity  dynamics,  scalar  mixing  and  combustion  in  high 
Reynolds  number,  high  Damkohler  number  flows  using  the  transport  element  method,  a 
Lagrangian  approach  to  the  simulation  of  the  Navier-Stokes,  species  and  energy 
conservation  equations,  has  several  advantages,  including  the  reduction  of  the 
computational  effort  by  concentrating  the  computational  elements  at  zones  of  high  vorticity, 
the  direct  evaluation  of  the  correct  (non  diffuse)  material  interface  across  which  reaction 
occurs,  and  the  natural  temporal  and  spatial  adaptivity.  Results  of  a  large  number  of 
simulations  which  we  have  conducted  over  the  past  few  years  have  shown  that: 

(1)  Mixing  and  combustion  in  a  shear  layer  depends  strongly  on  the  inlet  conditions  since 
a  small  shift  in  the  relative  initial  location  of  the  shear  (vorticity)  layer  with  respect  to 
the  mixing  (material  interface)  layer  can  lead  to  substantial  changes  in  the  striation  and 
whorling  of  the  reaction  surface,  as  shown  in  figure  1.  This  leads  to  substantial 
changes  in  the  amount  of  mixed  fluid  and  rate  of  product  formation  in  the  layer,  as 
shown  in  figure  2.  The  results  indicate  while  one  may  be  able  to  control  the  mixing 
rate  by  adjusting  judiciously  the  boundary  conditions,  accurate  predictions  require 
ensemble  averaging,  even  under  ergodic  conditions,  and  not  temporal  or  spatial 
averaging  (see  Pubs.  [7,9]). 

(2)  The  impact  of  the  boundary  conditions  on  the  mixing,  and  hence,  reaction  rates  extend 
to  the  case  when  heat  release  is  substantial.  This  is  because  when  the  reaction  smface 
is  located  further  away  from  the  vorticity  layer,  the  impact  of  heat  release  on  the  latter 
is  much  more  moderate  than  in  the  case  when  the  two  layers  are  aligned  This  issue 
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is  discussed  in  detail  in  the  paper  in  Appendix  I  where  we  investigated  means  by 
which  the  negative  impact  of  heat  release  on  the  mixing  rate  can  be  mitigated  by  (a) 
adding  a  wake  component  to  the  initial  shear  layer  profile  (boundary  layer  tripping)  or 
(b)  displacing  the  two  layers  with  respect  to  each  other. 

(3)  Even  at  relatively  low  Reynolds  numbers,  the  flame  zone  is  much  thinner  than  the 

large  scale  vorticity  stmcture  (characterized  by  AU  and  A,  the  shear  and  initial 
instability  wave  length,  respectively,  and  thinner  than  the  small  scale  foldings  inside 
these  stmctures  (described  by  the  striation  a =AU/e  where  eis  the  prevailing  strain). 

The  two-way  cascade  (formation  of  large  scale  due  to  rollup  and  pairing,  and  small 
scales  due  to  stretching  and  spiraling)  leads  to  the  formation  of  vortex  worms,  while 
the  flame  surface  continues  to  fold  inside  and  around  these  stmctures. 

(4)  The  flame  stmcture  within  these  folds  is  strongly  affected  by  the  competition  between 
convection  (stretch)  and  diffusion  normal  to  the  flame  surface.  Locally,  the  stmcture 
resembles  that  of  a  one-dimensional  stretched  flame,  except  when  local  flame  surface 
ciu^ature  is  much  higher  than  the  flame  thickness,  a  condition  encountered  rarely.  In 
cases  when  the  initial  material  interface  does  not  coincide  with  the  vorticity  layer, 
when  a  strong  wake  component  is  present  in  the  shear  layer  profile,  or  when  the  heat 
release  rate  is  rather  high,  the  flame  surface  may  exist  outside  the  large  vortices. 

(5)  The  reaction  zone  embedded  within  the  flame  stmcture  is  much  thinner  than  the 
convection-diffusion  zone  (the  ratio  of  the  two  is  the  Zeldovich  number).  Resolving 
the  reaction  zone  is  rather  demanding,  especially  at  strains  closer  to  the  extinction 
strain.  The  difference  between  the  two  thicknesses  accounts  for  the  large  number  of 
computational  element  required  and  the  enormous  reduction  of  the  time  step  needed 
during  finite  rate  reaction  simulations. 

(6)  While  along  most  of  the  flame  surface,  the  flame  acts  as  a  free  standing  flame,  i.e., 
located  between  two  free  streams  existing  far  away  of  the  reaction  zone,  within  the 
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folds  and  near  the  spifal  tip  opposite  flames  may  interact  in  such  a  way  that  their  local 
structures  differ  from  that  of  a  free  flame. 

1.4.  Implications  and  the  Elemental  flame  Model 

The  implications  of  the  above  observations,  supported  by  the  experimental 
observation  on  the  winding  of  material  interfaces  in  shear  flows  shown  in  figtires  3  and  4, 
include: 

(1)  Chemical  kinetics,  at  a  reduced  level,  is  necessary  in  modeling  the  chemical  structure 
(reaction-diffusion  zone)  of  the  flame,  and  especially  for  highly  strained  elemental 
flames. 

(2)  Beyond  their  birth  (introduction  at  the  zone  where  the  original  interface  stretches  the 
most),  elemental  strained  flames  continue  to  bum  as  they  move  along  the  spiral.  Their 
stmcture  changes  via  the  action  of  the  unsteady  strain  and  more  significantly  due  to  the 
changes  in  the  boundary  conditions  (within  folds  and  near  the  spiral  trip). 

(3)  At  the  spiral  epicenter,  elemental  flames  die  as  the  strain  vanishes  and  one  reactant  gets 
fully  consumed.  What  remains  inside  the  main  stmcture  at  the  eddy  center  is  a  mixture 
of  products  and  one  of  original  reactants  (depending  on  the  original  concentrations  of 
reactants  in  the  free  streams  and  the  stoichiometry)  which  is  not  capable  of 
communicating  with  the  free  streams.  (This  can  be  considered  as  the  initial  state  of 
formation/destraction  of  minor  or  secondary  species.) 

This  new  insight  led  us  to  the  formulation  of  the  elemental  flame  concept  as  a 
physically  acceptable  and  direct-simulation  supported  model  of  reaction  zone  stmcture  in 
turbulent  combustion  at  high  Damkohler  number.  The  concept  is  illustrated  graphically  in 
figures  5  and  6  where  both  zero  order  (zero  thickness  flame)  and  first  order  models  of 
premixed  and  diffusion  flames  are  shown.  In  this  model,  elemental  flames  are  introduced 
continuously  where  the  stain  rate  is  non-zero.  As  the  elemental  flames  move  closer  to  the 
large  vortical  stmcture,  and  as  their  stmctures  are  strained,  their  thickness  reach  an 
asymptotic  limit  of  Vwe.  The  spiraling  of  the  flames  around  the  stmcture  and  the 
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concomitant  stretch  of  its  structure  helps  maintain  its  vigor  by  removing  the  products  and 
bringing  in  reactants  closer  to  the  reaction  zone.  Finally  the  elemental  flame  joins  the 
central  region  of  the  eddy  where  it  contributes  to  the  accumulation  of  products  in  the 

macroelement  whose  radius  is  close  to  A.  Only  secondary  species  chemistry  can  be 

sustained  within  this  macroelement  element  The  results  of  the  implementation  of  this 
model  are  shown  in  figures  7  and  8  where  a  comparison  between  a  direct  simulation  (top), 
an  elemental  flame  model  at  high  resolution  (middle),  and  elemental  flame  model  at  low 
resolution  (bottom)  is  shown.  Clearly,  the  concept  is  capable  of  reproducing  the  results  of 
the  direct  simulation  faithfully. 

1.5.  Large  Stmcture  Simulations 

Based  on  direct  simulations  of  shear  layers  using  vortex  methods,  and  extending 
the  concepts  developed  to  treat  mixing  and  reaction  in  elemental  flames,  we  have 
formulated  a  new  algorithm  to  improve  the  efficiency  of  vortex  simulations  of  reacting 
flows,  without  sacrificing  any  of  its  inherent  advantages.  In  this  algorithm,  computational 
vortex  elements  represent  elemental  vortices  introduced  at  zones  of  strong  vorticity,  i.e., 
non-zero  shear,  whose  scales  are  representatives  of  the  local  vorticity  thickness.  Elemental 
vortices  are  then  strained  and  elongated  via  the  prevailing  strain  field  forming  vorticity 
microelements  as  they  spin  around  the  spiral  centers.  Scales  normal  to  the  vorticity  are 
determined  by  striation  due  to  folding.  Finally  these  micro  elements  are  combined  at  the 
spiral  center  to  form  vorticity  macro  elements.  These  macro  structures  continue  to  grow  by 
absorbing  micro  element  as  they  approach  their  centers,  by  diffusion,  break  up  through 
filamentation,  and  stretch  in  the  vorticity  direction.  Equilibrium  is  soon  followed  by 
irrelevance  (death)  as  the  vortical  structures  diffuse.  By  construction,  as  the  resolution  is 
refined,  this  algorithm  converges  to  the  original  vortex  method.  The  implementation  of  this 
concept  which  we  dubbed  as  large  structure  simulation  is  illustrated  in  figures  9  and  10 
where  the  top  part  of  each  figure  shows  results  of  direct  simulations,  the  middle  part  shows 
large  structure  simulation  with  high  resolution  and  the  bottom  part  depicts  large  structure 
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simulation  at  low  resolution.  The  figures  show  that  the  recovered  vorticity  field  in  all  three 
cases  is  nearly  the  same  while  the  number  of  elements  used  has  been  reduced  by  a  factor  of 
two,  reducing  the  computational  effort  by  a  factor  of  four. 

1.6.  Elemental  Flames  As  Reaction  Zone  Models  at  Infinite  Reaction  Rates 

The  implementation  of  the  elemental  flame  model  in  this  large  structure  simulation 
of  the  vorticity  field  to  the  case  of  infinite  chemical  kinetics  is  essentially  the  same  as 
described  before,  i.e.,  introduction  of  elemental  flames  at  the  interface  between  the  reacting 
species,  the  straining  and  elongation  of  these  elemental  flames  as  they  spiral  around  the 
center  of  the  macro  vortex  elements  leading  to  their  multiplication,  the  death  of  the 
elemental  flames  as  they  reach  the  tip  of  the  spirals  at  or  near  the  centers  of  the  vortex 
macro  element  forming  combustion  macro  element.  The  accumulation  of  burnt  mixture  at 
the  spiral  center  leads  to  the  formation  of  a  mixing  macro  element  which  resembles  a 
radially  stratified  reactor  whose  limit  of  perfect  homogeneity  is  a  well  stirred  reactor.  We 
note  that  using  a  Shvab-Zeldovich  formulation  of  the  combustion  problem  is  a  natural 
choice  here  since  the  location  of  the  elemental  flames  can  always  be  found  at  a  fixed,  e.g., 
zero  value  of  the  conserved  scalar  concentration.  This  is  shown  schematically  in  figure  12. 

1.7.  Elemental  Flames  as  Reaction  Zoner  Models  at  Finite  Rate  Reactions 

To  extend  the  above  concepts  to  the  case  of  finite  rate  kinetics,  we  combine  the 
elemental  flame  concept  with  the  solution  of  the  unsteady  strained  flame  described  in  detail 
in  Appendix  n  and  Appendix  HI.  In  both  papers  we  describe  a  model  in  which  the  flame 
structure  is  computed  under  the  assumption  that  the  stretch  rate  within  the  flame  is  uniform 
thus  by  reducing  the  complexity  of  the  problem  enormously.  In  this  case,  the  stracture  of 
elemental  flame,  in  terms  of  the  species  and  temperature  profiles,  under  the  action  of  the 
strain  exerted  by  the  layer,  shown  in  figure  1 1,  is  computed  independently  and  used  in  the 
immediate  vicinity  of  that  particular  elemental  flame  in  the  physical  domain  to  reconstruct 
the  overall  reaction  field.  The  luiion  of  these  elemental  flames  in  the  physical  domain 
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constitutes  the  reaction  zone  structure  in  the  turbulent  field.  The  implementation  of  this 
concept  in  the  shear  layer  is  shown  in  figure  13. 

Results  shown  in  figures  13  and  14  represent  the  first  implementation  of  the 
concept  of  finite  rate  kinetics  elemental  flames  in  a  turbulent  combustion  simulations.  As 
mentioned  above,  to  perform  these  simulations,  elemental  flames  are  introduced  wherever 
an  interface  between  reacting  species  exists  in  the  domain  (birth  of  an  elemental  flame). 
These  elemental  flames  are  strained  and  elongated  by  the  underlying  shear  thus  by 
multiplying  and  forming  an  ever  increasing  set  of  unsteady  strained  elemental  flames 
(USEF)  around  the  spiraling  areas  of  the  vortex  macroelements.  The  internal  structure  of 
the  elemental  flames  are  computed  by  integrating  the  uniformly  strained  flame  structure 
equations  as  in  Appendices  II  and  HI  using  a  single  step,  a  reduced  multistep  or  a  detailed 
kinetic  mechanism.  The  structure  calculations  are  then  mapped  back  over  to  the  physical 
domain  to  construct  the  reacting  field  there.  Straining  can  lead  to  local  quenching. 

At  the  tip  of  the  spiral  inside  the  macro  structure,  the  elemental  flames  are  starved  of 
one  of  the  primary  reactants,  thus  by  extinguished  (death).  Meanwhile,  elemental  flame  are 
spiraling  inwards  towards  the  center,  they  communicate  with  neighboring  elements  via  the 
boundary  conditions,  a  situation  which  could  lead  to  strong  changes  in  their  structure. 
Figure  14  shows  the  reacting  field,  represented  by  the  product  distribution,  of  a  shear  layer 
as  computed  using  this  concept  (top)  and  using  direct  simulations  (bottom).  Results  show 
that  the  model  reproduces  the  direct  simulation  results  very  well. 

1.8.  Uniform  Strain  Flame  Stmcture  Model  for  the  Reaction  2^ne 

The  uniform  strain  flame  model,  described  in  detail  in  Appendices  n  and  IE,  is 
based  on  the  assumption  that  there  exists  a  value  of  the  strain,  which  we  call  the  effective 
strain  which,  assumed  to  act  uniformly  throughout  the  flame  stmcture,  will  enable  us  to 
predict,  accurately,  the  species  and  temperature  profiles,  and  burning  rate,  without  having 
to  solve  the  momentum  equation.  We  have  found  that  using  this  approximation,  one  can 
transform  the  original  equations  into  a  set  of  reaction-diffusion  equations,  i.e.,  remove  the 
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troublesome  convective  derivatives,  and  thus  simplify  the  numerical  integration 
enormously.  The  equivalent  effective  strain  has  been  shown  to  be  the  external  applied 
strain  adjusted  to  correspond  to  the  reaction  diffusion  zone,  i.e.,  to  the  zone  where  the  . 
temperature  is  close  to  the  flame  temperature.  Using  this  model,  we  have  shown  that 
results  obtained  using  this  model  agree  very  well  with  the  solution  of  the  exact  equations, 
over  a  wide  range  of  strain  and  for  both  premixed  and  diffusion  flames.  We  have  also 
implemented  a  four  step  reduced  chemical  kinetic  scheme  for  methane  in  the  computation 
and  showed  that  it  works  well  for  both  the  premixed  and  diffusion  flame  cases. 

1.9.  Conclusions  and  future  Work 

The  work  summarized  in  this  report  describes  the  formulation  and  implementation 
of  an  efficient  and  accurate  reaction  zone  model  of  combustion  in  turbulent  flows  at  high 
Reynolds  and  Damkohler  numbers.  The  model  has  been  developed  to  mimic  the  observed 
structure  of  reaction  zones  in  turbulent  flows  obtained  by  Lagrangian  direct  simulations 
under  the  same  conditions.  In  particular,  the  reaction  zone  is  considered  locally  as  an 
unsteady  strained  flat  flame  whose  structure  is  obtained  by  integrating  the  stretched  flame 
equations  assuming  that  the  strain  rate  remains  constant  throughout  its  thickness.  The 
flame  is  affected  locally  by  the  folding  and  spiraling  induced  by  the  vorticity  structures 
since  the  boundary  conditions  used  in  integrating  its  equations  are  adjusted  to  reflect  the 
proximity  of  the  “free”  streams  to  the  reaction-diffusion  zone.  Moreover,  the  constant 
strain  is  taken  to  be  the  effective  strain  acting  at  the  same  zone,  to  improve  the  accuracy  of 
the  approximating.  Properties  of  folding  and  spiraling  material  elements  are  exploited  in 
removing  elemental  flames  as  they  move  very  close  to  each  other,  or  when  they  reach  the 
tip  of  the  spirals  at  the  center  of  the  large  vortical  stmcture.  Solutions  obtained  using  this 
model  compare  favorably  with  direct  simulations. 

Several  important  fundamental  issues  remain  to  be  addressed.  The  reaction  zone 
model  developed  in  this  work  applies  accurately  at  the  high  Reynolds  number  and 
Damkohler  number  limit.  As  we  have  seen  from  the  results  of  the  direct  numerical 
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simulation  results,  even  at  this  limit,  local  areas  of  nearly  homogeneous,  or  radially 
stratified  reactors,  could  exist  within  the  cores  of  the  large  vortical  stmctures  due  to  the 
rapid  accumulation  of  combustion  products  and  the  reduction  of  the  strain  rate  there.  The 
role  of  this  these  reactors  is  most  essential  in  determining  the  post  flame 
formation/destruction  of  minor  species,  such  as  CO  and  Nox.  Work  on  modeling  the 
dynamics  of  these  reactors  using,  e.g.,  well  or  partially  stirred  reactors,  would  generalize 
the  application  of  the  current  model  to  a  wider  range  of  governing  parameters. 

Another  significant  issue  the  behavior  of  the  elemental  flames  in  a  three  dimension 
space.  The  two  most  significant  problems  in  three  dimensions  are;  (1)  Developing  an 
efficient  approach  to  compute  the  evolution  of  a  material  interface  at  scales  comparable  with 
the  diffusive  scales.  This  problem  here  is  complicated  by  the  extreme  convolution 
experienced  by  material  elements  in  three  dimensions.  (2)  The  possibility  of  encountering 
interacting,  non  parallel  elemental  flames.  In  this  case,  it  is  not  obvious  when  one  can 
consider  neighboring  elemental  flames  as  being  independent  and  when  transport  fluxes 
along  the  flames  are  negligible.  To  maintain  computational  efficiency,  one  must  carefully 
eliminate  elemental  flames  as  they  “collide”  forming  locally  mixed  reacting  zones,  and 
develop  reaction  zone  models  for  these  zones.  (We  note  here  that  in  two  dimensions,  the 
probability  of  formation  of  such  zones  is  much  lower  since  flames  can  not  collide,  they  can 
only  meet  “face-to-face”  as  described  earlier  in  the  report.) 
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maximum  vorticity  and  maximum  scalar-gradient  on  tlie 
evolution  of  the  downstream  product  mass-fraction 
field.  The  misalignment  increases  from  top  to  bottom 
as  Od,  Id,  2d,  where  d=0.0585. 


9661  ^ZS6‘6P6  dd  ‘9  *££  *f  WIV  ‘^IOOOUbh  Q  H  H  X* 


**ZH  03  vs  tioAiJp  joXbi 

jBoqs  ‘[Buoisudimp-OMJ  b  jo  saS^un  po^tooi  osbiu  •£  ojtiSiii 


33.3  fns  41.7  ms 


Figure  4.  Definition  of  interfaces  and  areas  from  contour  analysis. 
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maximum  tip  angle  a=23  (middle)  and  a=47  (bottom). 


U  ' — ■ 

o  w  ^  S 

^  El  ^ 

S 

®  ® 

«£  S- 

-.=  JM 

S  ®  — 

w 

2  ®  JS  s 

f  tt « 

S  »£« 

aS  S  j£ 

CZ  PBi  W 

^  c<  ^  S 

05  s  - 

s^Jin 

•5  g  P  II 
w  «  «  rt 

^  o 

O  s  W 

</^  S  _  s 

v/5  ^  S  {5 

«  —  © 

SW  •!>«  s^ 

©  ^  .s^ 

<:■  'W 

■w  ■■■  *■> 

w  S  SlJ 

s  ^  -2  " 

*s  S  S  2 

2^2  = 
s-gj-  g 
^  S  ®  s 


00 

(U 

v-i 

b£) 

•  rH 


TIME- 10.00  STEP-*..j  ELEMENTS-  2493 


JS  iC 


.s  II 

2  T3 

•S  c 
S  o 

C3 

^  . 

O  C2  ^ 
^  C  C5 

WW  o 

"a  1-H  ^ 
one 
js  o  — • 

0$  .52 

S  CO 

(D  ^ 

•4-* 

c  S  3> 

Ct3  C 

St-«  K 

c3  ji 

^  ex,  *5 

«  j-  o 

X  >  -S 

Si 

^  c  o 
O  to  u 

>  o  B 

(D  x:  G 

^  ^  p. 

C4)  O  O 

C  60  *5 

'G(D 
CO  S'  *-1 

G  ^  G 

I  >9 

ii« 
^  r  § 

o  o 

»<*::  2 
5^  i2  « 
is  G  jG 
o  <u  ■*- 

>  S  o 

«t>  «  4-. 

J3  ^  t) 
■“  fci  4> 
«4-,  o  a 

o  2  «? 

c  £  ^ 
•2  >^’§ 
3o  > 

»f-H 

>  «- 

«  5  2 

^  ^  CJ 
KJ  ^  « 
•iS  o  > 

C^  4-»  >. 

Cl,  Ch^ 
oo  <D  *3 

CJ  X 

s  ® 

H  o  > 


inetluKi  (top),  aiici  that  incorporating  the  concept 
of  >orticity  macroelements  for  two  values  of  the 
growth  parameter  c=l  (middle)  and  c=2  (bottom). 


of  the  material  line  shown  above 


The  infinite  reaction  rate  approximation  of  the  diffusion  flame; 


Product  mass-fraction  obtained  from  the  flamelet 
solutions  for  non-iiityeracting  flamelets  (top)  and 
interacting  flamelets  (bottom). 
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ABSTRACT 

The  effects  of  the  presence  of  a  wake  component 
in  the  inlet  velocity  profile,  on  the  dynamics  of  a  non- 
premixed  exothermic  shear  layer  are  investigated 
numerically.  The  Lagrangian  transport  element  method  is 
used  to  obtain  two-dimensional  simulations  of  the  flow. 
Combustion  is  modeled  using  an  infinite  reaction  rate 
model  based  on  the  Shvab-Zeldovich  formulation. 
Simulations  with  and  without  combustion  exothermicity 
are  obtained  for  two  different  inlet  conditions:  One  with 
an  errorfunction  inlet  velocity  profile,  and  one  in  which 
this  profile  includes  a  Gaussian  wake  deficit.  Results 
indicate  that  the  presence  of  the  wake  component  in  the 
inlet  velocity  profile  alters  the  downstream  flow  and  its 
interaction  with  the  exothermic  reacting  field  due  to  two 
main  reasons:  (i)  The  wake-modified  profile  is  much  more 
unstable  than  the  errorfunction.  As  a  result  it  yields 
higher  downstream  rates  of  mixing  and  burning,  (ii)  The 
relative  location  of  the  reaction  interface  with  respect  to 
the  vorticity  field  differs  between  the  errorfunction  and 
the  wake-modified  cases.  As  a  result,  even  though 
combustion  exothermicity  decreases  the  burning  in  both 
cases  as  compared  to  their  corresponding  non-exothermic 
flow,  it  does  so  via  different  mechanisms.  In  the 
errorfunction  flow  the  reaction  interface  is  located  at  the 
center  of  the  single  sign  vorticity  layer  which 
characterizes  the  vorticity  field.  The  reduction  in  the 
burning  in  this  case  is  a  result  of  the  decay  of  the  vorticity 
field  due  to  volumetric  expansion  which  delays  the  rollup 
and  results  in  weakened  downstream  vortical  structures. 
In  the  wake-modified  case  the  reaction  interface  is  located 
between  two  vorticity  layers  of  opposite  signs.  The  layer 
with  the  same  sign  as  that  of  the  errorfunction  case 
dominates  the  dynamics.  In  this  case,  however,  the 
reaction  interface  is  initially  located  outside  it  and  gets 
wrapped  around  it  as  the  this  layer  gets  destabilized. 
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The  reduction  in  the  burning  in  the  exothermic  case 
comes  about  via  a  delay  in  this  interaction  of  the  vorticity 
layer  and  the  reaction  interface  as  volumetric  expansion 
pushes  the  two  away  from  one  another.  The  same 
mechanism  is  also  responsible  for  a  delay  in  the 
interaction  of  the  two  vorticity  layers  of  the  wake- 
modified  flow  with  the  consequence  that  the  wake  region 
survives  further  downstream. 


I.  INTRODUCTION 

Numerical  simulations  of  turbulent  reacting 
shear  layers  in  which  the  governing  equations  are 
integrated  directly  without  averaging  and/or  resorting  to 
turbulence  closure  modeling  have  been  commonly  used  to 
elucidate  the  mechanisms  of  interaction  between 
turbulence  and  combustion  in  shear  flows  [e.g.  1-6]. 
These  flows  are  dominated  by  large  scale  vortical 
structures  [7,  8]  whose  presence  and  impact  on  the  mixing 
and  combustion  cannot  be  directly  captured  using 
turbulence  closure  and  existing  turbulence-combustion 
interaction  models,  respectively,  and  hence  must  be 
directly  computed. 

Recently,  these  simulations  have  enabled  the 
systematic  and  detailed  investigation  of  the  effect  of 
combustion  exothermicity  on  the  flow  physics  [4-6].  In 
this  regard,  it  has  been  shown  that  the  mixing 
enhancement  by  the  enlargement  of  the  flame  surface  area 
due  to  the. action  of  the  large  scales  plays  a  key  role  in 
burning  augmentation.  Moreover,  it  has  also  been 
concluded  that  exothermic  energy  impacts  the  dynamics 
in  such  a  way  that  the  overall  growth  rate  of  the  mixing 
zone  is  reduced  with  respect  to  that  observed  in  the  non- 
exothermic  flow.  This  reduction,  which  takes  place 
despite  the  increase  in  volume  due  to  the  deposition  of 
exothermic  energy,  is,  for  the  most  part,  a  result  of  the 
decay  of  the  vorticity  field  due  to  volumetric  expansion. 
It  has  also  been  established  that  the  structure  and 
interactions  of  the  vortices  are  altered  in  the  exothermic 
flow  due  to  the  redistribution  of  vorticity  by  the  baroclinic 
mechanism. 

The  validity  of  these  conclusions,  however,  is 
restricted  by  the  limiting  assumptions  used  to  obtain  the 
corresponding  results.  One  common  source  of  contention 
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is  the  simplification  of  the  inlet  conditions.  This  is 
particularly  troublesome  when  one  considers  the  well 
known  sensitivity  of  the  flow  to  these  conditions  [8,9].  In 
the  studies  noted  above,  the  inlet  velocity  profile,  which 
due  to  the  presence  of  the  upstream  splitter  plate 
separating  the  two  streams  includes  a  wake  component,  is 
smoothened  into  a  monotonically  varying  one,  usually  an 
errorfunction  or  a  hyperbolic  tangent.  The  origin  of  this 
assumption,  which  is  motivated  by  numerical 
considerations,  lies  in  the  study  of  the  non-reacting  flow. 
Early  experimental  studies  indicated  that  in  a  non-reacting 
flow  the  wake  component  of  the  mean  velocity  profile  is 
quickly  annihilated  downstream  of  the  splitter  plate  [7, 9]. 

Linear  stability  analysis  of  the  flow  also  points  in 
die  same  direction.  In  a  number  of  studies  it  has  been 
shown  that  in  contrast  to  a  single  mode  of  instability  of 
the  monotonically  varying  profile,  a  velocity  profile  with 
a  wake  component  exhibits  two  modes  of  instability, 
namely  the  ‘sinus’  or  ‘shear  layer’  and  the  ‘varicose’  or 
‘wake’  modes  [10-13].  In  the  vast  majority  of  conditions, 
however,  the  sinus  mode,  which  is  similar  to  that  of  the 
monotonic  profile,  albeit  with  a  different  growth  rate  and 
phase  speed,  was  found  to  dominate.  As  a  result,  one  may 
hope  to  capture  the  essential  features  of  the  flow,  at  least 
in  a  qualitative  sense,  by  using  the  monotonic  profile. 

Numerical  simulations  of  non-reacting  shear 
layers  using  a  monotonically  varying  inlet  profile  were 
successful  in  reproducing  the  flow  instantaneous  and 
mean  dynamics  [e.g  14,  15].  We  believe,  however,  that 
the  effectiveness  of  this  type  of  inlet  condition  in  the  non- 
reacting  flow  does  not  necessarily  imply  its  applicability 
to  the  exothermic  reacting  flow.  A  number  of  important 
issues  suggest  that  such  an  extension  may  be  unjustified. 
First,  in  the  reacting  case  a  variable  density  field  is 
present.  The  studies  of  the  linear  stability  of  the  flow 
mentioned  earlier  clearly  indicate  that  the  flow 
destabilization  is  modified  in  the  presence  of  a  variable 
density  field.  The  probability  of  exciting  the  varicose 
mode  may  increase,  particularly  when  the  thickness  of  the 
density  profile  is  much  smaller  than  that  of  the  velocity 
profile.  More  important,  however,  is  the  issue  of 
volumetric  expansion.  The  flow  arrangement  at  the  inlet 
implies  that  deposition  of  exothermic  energy,  and  hence 
volumetric  expansion,  will  take  place  inside  the  wake 
region.  The  possibility  that  this  may  have  substantial 
effects  on  the  downstream  dynamics  cannot  be 
disregarded. 

In  this  paper,  we  investigate  the  impact  of  the 
wake  component  of  the  inlet  velocity  profile  on  the 
exothermic  flow,  by  comparing  results  of  numerical 
simulations  in  which  such  an  inlet  profile  is  used,  with 
those  of  corresponding  simulations  with  a  monotonically 
varying  inlet  profile.  The  transport  element  method  is 
used  to  provide  high  resolution  two-dimensional 
simulations  of  the  flow.  The  reacting  field  is  modeled 
using  an  infinite  reaction  rate  model  based  on  Shvab- 
Zeldovich  conserved  scalars.  In  earlier  work  we  have 
shown  this  model  to  be  effective  in  enabling  the 
prediction  of  the  effects  of  combustion  exothermicity  on 
the  reacting  field,  as  long  as  the  reaction  is  fast  connpared 


to  the  flow  time  scale  [16].  Thus,  the  results  shown 
herein  ^ply  in  this  limit. 


n.  FORMULATION  AND  NUMERICAL  SCHEME 

The  details  of  the  formulation  and  numerical 
scheme  used  in  the  simulation  of  infinite  reaction  rate 
exothermic  shear  flow  are  given  in  Refs.[16, 17].  Herein 
only  a  summary  of  the  main  features  is  provided. 

Formulation 

A  two-dimensional  high  Reynolds  number,  non- 
premixed,  exothermic  reacting  shear  flow  is  considered. 
Compressibility  effects  are  allowed  only  in  the  low  Mach 
number  limit,  the  effects  of  gravity  are  neglected,  all 
species  are  assumed  to  behave  as  perfect  gases  with  equal 
molar  masses  and  constant  and  equal  diffusivities  and 
specific  heats,  and  combustion  to  occur  according  to  a 
single  step,  irreversible,  mole  preserving,  infinite  rate 
reaction. 

Under  the  infinite  reaction  rate  model  reactants 
bum  completely  upon  contact.  As  a  result,  the  reaction 
zone  collapses  onto  an  interface  which  separates  the  two 
reactants.  The  solution  of  such  a  reacting  field  which 
experiences  discontinuities  in  the  gradients  of  the 
primitive  scalar  variables,  is  facilitated  by  the  introduction 
of  Shvab-Zeldovich  (SZ)  conserved  scalars  which  are 
piecewise  continuous  in  the  whole  domain. 

The  resulting  governing  equations  are 
transformed  into  a  vorticity — ^scalar-gradient  form.  This 
is  accomplished  by  decomposing  the  velocity  field  into 
expansion,  potential  and  vortical,  components  using  the 
Helmholtz  decomposition  and  representing  the  first  two 
by  the  appropriate  form  of  the  continuity  equation  and  the 
third  by  the  streamfunction  equation.  The  latter 
necessitates  the  solution  of  the  vorticity  equation  which  is 
established  by  taking  the  curl  of  the  momentum  equation. 
For  the  scalar  field  the  transport  equations  of  the  gradients 
of  the  SZ  variables  are  established,  i.e.  in  non- 
dimesionalized  form. 


B  =  V4>,  +  V0p  +  Vj<  'm  ; 
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A  =  ri-0}'2  and  y=T-  ^Yp  ■,  (7a, b) 

1+0 

^=:.gVu-gx(  ak)+  -^V^g  ,  (8) 

dt  Pe 

where  g  =  VA  or  Vy  . 

In  the  above,  u  =  (u,v)  is  the  velocity  vector  in  an 
X  =  (jc,y)  right-handed  Cartesian  coordinate  system,  t  is 

time,  V  =  (-^,-^)  and  -^  =  — +«-V.  0e  is  th® 
\dxdyl  dt  dt 

expansion  velocity  potential,  <Pp  is  the  velocity  potential 
related  to  the  boundary  conditions,  ^  is  the 

streamfimction,  (ok  =  Vatu  is  the  vorticity,  where  k  is  the 
unit  vector  normal  to  the  plane  of  motion,  p  is  the 
pressure,  p  is  the  density  and  T  is  the  temperature.  Re  and 
Pe  are  the  Reynolds  and  Peclet  numbers,  respectively, 
-non-dimensionalizing  scales  are  given  in  Section  ni.  Yi 
is  the  mass-fraction  of  species  rii,  /  =  1, 2  for  the  reacting 
species,  i  =  p  for  the  product,  0*  and  0  are  the  molar  and 
mass  stoichiometric  ratios,  respectively,  and  Qo  is  the 
normalized  enthalpy  of  reaction. 

Finally,  it  is  clarified  that  the  reconstruction  of 
the  reacting  scalar  field  from  the  SZ  solutions  is 
accomplished  using  equations  (6)  and  (7)  and  the 
knowledge  that  reactants  cannot  coexist: 


ASO 

II 

F2  =  0,  Yp 

=  1-A, 

and 

T—  Y+  ( 

1+0 

1-A) 

(9a) 

A<0 

o 

II 

F2  =  - 

+ 

II 

A 

0’ 

and 

T—  Y+  ( 

1+0 

l  +  i^) 

0 

(9b) 

It  is  noted  that  the  above  equations  define  A  =  0  as  the 
locus  of  the  reaction  interface. 

Numerical  Scheme 

The  governing  equations  are  solved  using  the 
transport  element  method.  In  this  Lagrangian,  grid-free 
numerical  scheme  the  vorticity,  divergence  and  SZ 
gradient  fields  are  discretized  over  a  number  of  elements 
of  finite  area,  strength  and  overtyping  Gaussian  cores. 

Evolution  of  the  flow  and  scalar  fields  is 
achieved  by  advecting  the  elements  with  the  local  velocity 
vector  while  at  the  same  time  updating  their  properties  by 
numerically  integrating  the  governing  transport  equations. 


This  is  accomplished  in  two  fractional  steps.  In  the  first 
step,  which  includes  all  processes  other  than  diffusion,  the 
elements  are  transported  and  their  strength  is  adjusted  as  a 
result  of  the  numerical  integration  of  the  non-diffusive 
transport  equations.  For  the  vorticity,  the  use  of  the 
circulation  equation,  with  the  material  acceleration 
substituting  the  pressure  gradient  in  the  baroclinic  term,  is 
preferred  due  to  its  simplicity.  For  the  SZ  gradients,  the 
governing  equation  is  substantially  sinylified  using  flow 
kinematics  to  relate  the  gradient  evolution  to  that  of  the 
material  line  stretch.  The  latter  is  readily  available  due  to 
the  Lagrangian  nature  of  the  scheme. 

In  the  second  fractional  integration  step,  diffusion 
effects  are  simulated  using  the  core  expansion  scheme.  In 
the  current  implementation  of  this  scheme  every  step  of 
core  expansion  is  followed  by  a  rediscretization  step 
which  returns  the  cores  to  their  original  size  while 
modifying  the  element  strength.  This  avoids  problems 
associated  with  the  transport  of  large  cores. 

The  severe  distortion  of  the  flow-computation 
map  due  to  the  Lagrangian  transport  may  lead  to 
deterioration  of  the  accuracy  of  the  solution  due  to 
insufficient  core  overlap.  This  is  prevented  by 
continuously  inserting/removing  elements  in  regions  of 
high  tensile/compressive  strains,  as  assessed  by  the 
distance  between  neighboring  elements.  This 
insertion/removal  is  carried  out  according  to  local 
conservation  laws. 

The  velocity  field  at  any  given  time  is 
established  by  obtaining  the  vortical  and  expansion 
components  via  convolutions  over  the  fields  of  the  vortex 
and  expansion  elements,  respectively,  and  the  potential 
component  using  a  Schwartz-Christoffel  conformal 
mapping.  The  aforementioned  convolutions  result  from 
the  integral  solutions  of  the  corresponding  Poisson 
equations,  Eqs.(2a,2c),  in  an  unbounded  domain. 
Bounded  domains  are  simulated  using  the  concept  of 
images.  For  the  scalar  field  a  similar  convolution  is 
utilized  by  realizing  that  the  scalar  and  its  gradient  may 
also  be  related  by  a  Poisson  equation  (i..e.  =  V-g ). 


Slip  walls.  Elements 

impermeable  deleted 


Figure  1.  Schematic  representation  of  the  geometry  of  the 
computational  domain  together  with  the  initial 
element  configuration  and  some  of  the  boundary 
conditions 
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IV.  PROBLEM  STATEMENT 


A  twiMlimensional  channel  geometry  of  height 
H  and  length  Xmax  is  considered  as  shown  in  Fig.l.  Two 
parallel  fluid  streams  of  different  velocities  but  same 
density,  po,  and  temperature,  To,  previously  separated  by  a 
thin  splitter  plate  located  half  way  between  the  walls,  mix 
along  the  channel  length.  Each  stream  carries  a  single 
reactant.  H,  po  and  T„  are  used  as  the  reference  length, 
density  and  temperature,  respectively,  while  Uu  the 
velocity  of  the  top,  fast,  stream,  as  the  reference  velocity. 

The  top  and  bottom  walls  are  modeled  as  rigid, 
slip,  impermeable  and  adiabatic  planes.  At  the  exit, 
uniform  conditions  are  assumed  for  the  potential 
components  while  for  the  rotational-gradient  components 
the  corresponding  elements  are  simply  deleted  as  they  exit 
the  domain. 

Two  sets  of  inlet  conditions  are  considered, 
identiHed  as  the  ‘errorfunction’  and  the  ‘wake-modified’ 
conditions.  Their  main  difference  lies  in  that  the  first 
includes  the  effects  of  the  presence  of  the  wake  which  is 
reminiscent  of  the  boundary  layers  formed  along  the  walls 
of  the  upstream  splitter  plate,  while  the  second  does  not. 
This  is  in  accordance  widi  the  main  objective  of  this  p^r 
which  is  to  investigate  the  importance  of  the  presence  of 
the  wake  deficit  in  the  inlet  velocity  profile  on  the 
exothermic  flow  dynamics. 

The  two  sets  of  conditions  are  specified  under 
the  following  general  profiles  for  the  velocity  and  the  S-Z 
scalars: 


where  erf  and  exp  are  the  errorfunction  and  exponential 
functions,  respectively,  y*  =  y-  0.5,  (Ta  and  Og  are  the 
standard  deviations  of  the  related  vorticity  and  scalar 

gradient  profiles,  respectively,  r  =  ^  is  the  velocity  ratio, 

Ui 

W  is  the  wake  component  strength  and  y©  is  a  cross¬ 
stream  adjustment  length  which  allows  the  wake-modified 
profile  to  experience  its  minimum  at  y  *  =  0. 

For  both  conditions  r=:0.5  and  0=1.  For  the 
errorfunction  condition  both  W  and  y®  are  set  equal  to 
zero  while  for  the  wake-modified  condition  W=  1.1  and 

y„  =  ^  l^-^l .  The  thickness  of  the  velocity  profile  for 

W  \  2  I 

both  conditions  is  Om  =  Caw  -  =  0.039.  Both  these 

velocity  profiles  are  shown  in  Fig.2  which  displays  the 
profile  specified  by  Eq.(lO)  for  a  variety  of  values  of  W. 


Inlet  Velocity  Profile 


Figure  2.  The  velocity  profile  described  by  Eq.(lO)  for 
r  =  0.5,  Ca  -  0.039  for  various  values  of  the 
wake  parameter  W.  Cases  with  W  =  0  and  1 . 1 
are  used  as  inlet  conditions  in  the  numerical 
simulations. 


The  inlet  species  profiles 
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Figure  3.  The  species  mass-fraction  profiles  specified  at 
the  inlet  for  the  errorfunction  inlet  profile  flow. 


The  choice  of  the  SZ  scalar  profiles,  like  those  of 
the  velocity  profiles  is  motivated  by  the  flow  physics. 
Equations  (9)  imply  that  the  specification  of  the  A  inlet 
profile  as  an  errorfunction  suggests  that  the  reacting 
species,  in  the  regions  where  they  exist,  will  also  be 
errorfunctions.  As  a  result  the  product  follows  a  profile 
which  experiences  a  sharp  peak  in  the  neighborhood  of 
the  reaction  region.  The  choice  of  7  as  unity  implies  that 
temperature  variation  is  only  a  consequence  of  reaction. 
(It  is  noted  that  this  choice  together  with  the  nature  of  the 
governing  transport  equations  and  the  rest  of  the  boundary 
conditions  given  earlier  imply  that  =1.  As  a  result 
the  computational  solution  is  substantially  simplified.) 
These  scalar  profiles  which  are  shown  in  Fig.3  for  the 
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eiTorfunction  inlet  condition  for  which  the  thickness  of  the 
scalar  profiles  is  the  same  as  that  of  the  velocity  profiles, 
i.e.  OgE  -  0(0,  are  not  unlike  those  encountered 
downstream  of  the  splitter  plate.  For  the  wake-modified 
inlet  condition  the  scalar  profiles  are  assumed  to  be 
thinner  than  those  of  the  errorfunction  case,  and  hence 
than  the  velocity  profiles  as  well.  Specifically 
Ogv  =  0.5  0(0.  This  is  to  approximate  the  unmixed  nature 
of  the  scalar  profiles  just  downstream  of  die  splitter  plate. 
It  is  noted  that  the  particular  choice  of  this  thickness,  as 
well  as  the  earlier  choice  of  W,  are  limited  by  discretizing 
constrains.  The  errorfunction  scalar  profiles  are  thicker 
due  to  the  assumed  mixing  between  the  edge  of  the 
splitter  plate  and  the  location  of  the  inlet  condition  which, 
as  was  stated  earlier,  is  assumed  to  annihilate  the  wake 
conqionent  of  the  velocity  profile. 


(a) 


Growth  Rate  of  Linear  Waves 


(b) 


Phase  Speed  of  Linear  Waves 


Figure  4.  The  growth  rate  (a)  and  phase  speed  (b)  of 
spatially  growing  linear  instability  waves  of  a 
parallel  shear  flow  with  base  velocity  profile 
that  described  by  Eq.(lO).  Curves  for  different 
values  of  the  wake  parameter  W  are  indicated. 
The  particular  cases  shown  correspond  to  those 
ofFig.2. 


The  particular  choice  in  the  value  of  the 
thickness  of  the  inlet  velocity  profiles  noted  above  is 
motivated  by  the  stability  characteristics  of  the  flow. 
Equation  (10)  indicates  that  the  velocity  profile  chosen  is 
essentially  an  errorfunction  profile  with  a  Gaussian  wake 
deficit  conqionent.  Stability  analysis  of  a  flow  using  this 
profile  as  a  base  profile  indicates  that  the  wake 
component  alters  the  stability  characteristics  from  diat  of 
the  plain  errorfunction  profile  in  a  non-trivial  way.  First, 
it  introduces  a  second  instability  mode,  die  ‘varicose’  or 
‘wake’  mode,  to  the  single  mode  of  the  errorfunction,  i.e. 
the  ‘sinus’  or  ‘shear  layer’  mode  [10-13].  It  is  found, 
however,  that  under  most  conditions  the  sinus  mode  still 
dominates.  Second,  it  amplifies  the  growth  of  the  sinus 
mode.  This  is  shown  in  Fig.4  which  displays  the  spatial 
growth  rate  (a)  and  phase  speed  (b)  of  this  mode  as  a 
function  of  the  wake  strength,  W,  for  the  non-reacting 
uniform  density  flow.  The  cases  displayed  correspond  to 
the  velocity  profiles  of  Fig.2.  Clearly  the  magnitude  of 
the  growth-rate  of  the  instability  is  highly  sensitive  to  the 
wake  strength  W,  increasing  with  increasing  values  of  this 
parameter.  The  most  unstable  frequency  of  the  flow,  on 
the  other  hand,  does  not  seem  to  be  strongly  dependent  on 
W.  The  phase  speed  of  the  instability  waves  is  also 
strongly  modified  by  the  wake,  decreasing  with  increasing 
the  waice  deficit,  particularly  so  around  the  most  unstable 
frequency.  This  implies  that  the  wavelength  of  the  most 
unstable  mode  also  changes  significantly.  For  the  profiles 
considered  here,  i.e.  where  W  =  0.0  and  1.1  the 
corresponding  most  unstable  wavelengths  are 
Ae  =  12.8<Tib  and  Aw  =  6.5O(0.  That  is,  the  wavelength  of 
the  unstable  waves  arising  from  the  wake-modified  profile 
is  approximately  half  of  that  of  the  errorfunction  profile. 

Simulations  of  shear  layers  with  errorfunction 
inlet  profiles  indicate  that  the  wavelength  of  the  linear 
waves  remains  approximately  constant  as  the  flow 
develops  into  the  non-linear  stage,  i.e.  it  is  indicative  of 
the  size  of  the  first  generation  eddies  [15].  Thus,  one 
would  expect  the  eddies  of  the  wake-modified  profile  to 
be  about  half  the  size  of  those  of  the  errorfunction.  The 
choice  of  a  for  the  errorfunction  profile  noted  earlier  is 
such  that  two  such  eddies  fit  inside  the  channel  height. 
This  is  done  with  the  aim  that  a  paired  eddy  may  be 
accommodated  in  the  channel.  In  contrast,  for  the  wake- 
modified  case  this  implies  that  the  channel  may  fit  about 
four  eddies.  This  argument  also  implies  that  with  half  the 
domain  length  in  the  case  of  the  wake-modified  flow 
similar  number  of  eddies  may  be  captured.  For  this 
reason,  and  with  the  aim  of  limiting  the  computational 
cost,  in  the  simulations  the  length  of  the  domain  was 
varied  between  Xmax  =  8  and  4  for  the  errorfunction  and 
the  wake-modified  inlet  profiles,  respectively. 

Initialization  of  the  calculation  is  carried  out  by 
assuming  that  the  inlet  conditions  persist  throughout  the 
domain.  Hence,  transport  elements  are  distributed  over 
nine  flat  material  layers  (lines)  lying  within  the  region  of 
finite  vorticity  and  A-gradient  defined  by  the  above 
profiles.  The  elements  are  of  square  area  of  side 
A=0.0195.  The  value  of  the  core  radius  is  5=0.0234,  i.e. 
A/5=0.833.  The  time-step  in  all  simulations  was  4r  =  0.1 . 
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V.  RESULTS  AND  DISCUSSION 


errorfunction  case,  in  the  wake-modified  case  the  shear 
and  mixing  regions  do  not  always  coincide. 


The  main  objective  of  this  paper  lies  in  analyzing 
the  effects  of  the  presence  of  the  splitter  plate  related 
wake  component  in  the  mean  velocity  profile  on  the  flow 
dynamics  in  the  case  where  exothermic  combustion  is 
present.  To  achieve  this,  simulations  at  a  fixed  Reynolds 
and  Peclet  number.  Re  =  Pe  =  12800,  with  and  without 
combustion  exothermicity,  Qo  =  0  and  6,  were  carried  out 
for  both  the  errorfunction  and  the  wake-modified  inlet 
condition  flows. 

Instantaneous  visualizations  of  both  the 
errorfunction  and  the  wake-modified  simulation  results  in 
terms  of  the  vordcity  and  product  mass  fraction  are  shown 
in  Figs.  5  and  6,  respectively.  Parts  (a)  and  (b)  of  each 
figure  display  the  errorfunction  inlet  profile  flow.  In  part 
(a)  the  entire  computational  domain  is  shown  (i.e. 
Xmox  =  8)  while  in  part  (b)  a  domain  equal  in  size  to  that 
of  the  wake-modified  flow  is  presented  (i.e.  Xmax  =  4)  - 
note  that  the  size  of  this  domain  is  also  indicated  in  part 
(a).  Part  (c)  displays  the  wake-modified  inlet  condition 
flow  results.  Each  part  con:^)ares  the  Qo  =  0  (top)  with  the 
Qo  =  6  (bottom)  cases.  The  locus  of  the  reaction  interface 
is  indicated  in  all  visualizations  by  a  daiic  line. 

V.l  Overall  Features 

Figures  5  and  6  show  that  there  are  differences  as 
well  as  similarities  between  the  errorfunction  and  wake- 
modified  flows.  Two  major  conclusions:  (i)  The  wake- 
modified  flow  gets  destabilized  earlier  and  more 
vigorously  than  its  errorfunction  counterpart.  This  is  in 
agreement  with  the  stability  analysis  results  of  the 
previous  section.  As  a  result,  the  convolution  of  the 
reaction  interface  is  much  more  significant  for  the  wake- 
modified  flow  and  the  resulting  burning  is  clearly 
enhanced,  (ii)  For  both  the  errorfunction  and  wake- 
modified  flow,  exothermicity  tends  to  delay  the 
destabilization  of  the  flow  and  the  intensity  of  the  rollup. 

Both  conclusions  are  confirmed  by  the  mean 
flow  results  shown  in  Fig.7,  where  the  growth  of  the  mean 
shear  region  (a)  and  of  the  product  integral  (b),  for  both 
the  errorfunction  and  the  wake-modified  flows,  are 
depicted.  The  shear  region  is  based  on  the  streamwise 
evolution  of  the  mean  velocity  profile.  It  is  defined  as  the 
distance  between  the  points  of  the  profile  where  the 
velocity  first  deviates  by  0.01(1%  of  t/i)  from  the  value 
of  the  neighboring  free  stream.  The  product  integral  is 

defined  as  Ip  =|  Yp  dy.  It  is  clear  that  the  shear  region 

and  the  amount  of  product  formed  are  substantially 
diminished  when  the  errorfunction  inlet  profile  is  used.  It 
is  also  seen  that  for  both  the  errorfunction  and  wake- 
modified  flows,  exothermicity  reduces  the  amount  of 
product  formed  over  most  of  the  domain.  A  similar 
conclusion  cannot  be  drawn  for  the  mean  shear  region, 
however,  particularly  in  the  case  of  the  wake-modified 
flow  which  exhibits  a  more  complex  behavior.  As  will  be 
clarified  in  what  follows  this  is  because,  unlike  the 


(a) 


Mean  Shear  Region  Thickness 


Figure  7.  The  mean  shear  region  (a)  and  the  product 
integral  (b)  the  errorfunction  and  wake- 
modified  flows  with  and  without  combustion 
exothermicity. 


V.2  Detailed  features 

The  instantaneous  realizations  shown  in  Figs.  5 
and  6  can  be  used  to  answer  several  important  questions 
relating  to  the  mechanisms  by  which  the  reduction  in  the 


Figure  5.  (Figure  on  opposite  page)  Instantaneous 
vorticity  distributions  for  the  errorfunction  (a- 
b)  and  the  wake-modified  (c)  flows.  Part  (a) 
displays  the  entire  errorfunction  flow 
computational  domain  while  part  (b) 
concentrates  on  a  section  equal  in  size  to  the 
domain  of  the  wake-modified  flow.  Each  part 
of  the  figure  contrasts  the  j2o  =  0  (top)  and 
Qo  =  6  (bottom)  cases.  In  all  cases  the  locus  of 
the  reaction  interface  is  indicated  (black  line). 
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Figure  5.  Caption  on  opposite  page 


Figure  6.  Caption  on  following  page 


Figure  6.  (Figure  on  previous  page)  Instantaneous 
product  mass-fraction  distri-butions  for  the 
errorfunction  (a-b)  and  wake-modified  (c) 
flows.  The  distributions  depicted  correspond 
to  those  of  the  vorticity  of  Fig.5.  Part  (a) 
displays  the  whole  errorfunction  flow 
computational  domain  while  part  (b) 
concentrates  on  a  section  equal  in  size  to  the 
domain  of  the  wake-modified  flow.  Each  part 
of  the  figure  contrasts  the  =  0  (top)  and 
Go  =  6  (bottom)  cases.  In  all  cases  the  locus  of 
the  reaction  interface  is  indicated  (black  line). 


burning  due  to  combustion  exothermicity  comes  about  in 
the  two  flows.  For  instance,  can  one  claim  that  the  only 
difference  between  the  wake-modified  and  errorfunction 
cases  is  the  one  related  to  the  streamwise  location  and  the 
growth  rate  of  the  initial  destabilization?  And  if  that  is 
the  case,  can  the  errorfunction  flow  be  used  as  a  good 
approximation  of  the  wake-modified  behavior  via  a 
spatial  shifting  or  scaling?  The  figures  indicate  that  this  is 
clearly  not  true.  The  errorfunction  and  wake-modified 
flows  exhibit  differences  in  their  evolution  of  fundamental 
nature.  Examples  of  these  are  the  origin  and  locus  of  the 
positive  (conterclockwise  rotating)  vorticity  and  the 
relative  location  of  the  reaction  interface  with  respect  to 
the  vortical  structures. 

V.l.LThe  non~exothermic,  uniform  density  flow 
The  vorticity  field 

Even  in  the  absence  combustion  exothermicity 
and  of  the  associated  density  variation,  the  evolution  of 
the  wake-modified  flow  exhibits  important  differences 
from  the  errorfunction  flow.  When  the  density  is  uniform, 
the  vorticity  of  a  fluid  element  may  change  only  via 
diffusion.  As  a  result,  its  sign  caimot  change.  Thus,  the 
sign  of  the  vorticity  field  in  the  non-exothermic  flow  is 
uniquely  determined  by  the  inlet  condition.  For  the 
errorfunction  flow,  where  the  inlet  velocity  profile  is 
monotonic,  the  inlet  vorticity  is  of  a  single  sign,  namely 
negative  (top  stream  faster).  Hence  the  vorticity  of  the 
downstream  field  is  also  always  negative.  In  contrast,  for 
the  wake-modified  flow  the  non-monotonicity  of  the  inlet 
profile,  results  in  vorticity  of  both  signs.  These 
observations  are  evident  in  Fig.5  which  also  points  to  their 
significance  to  the  downstream  flow  dynamics. 

For  the  errorfunction  case,  the  evolution  of  the 
flow  is  manifested  through  the  destabilization  and 
convolution  of  a  single  sign  vorticity  layer  which  leads  to 
the  formation  of  a  series  of  vorticd  structures.  For  the 
wake-modified  case,  on  the  other  hand,  the  flow  evolves 
through  the  interaction  of  two  vorticity  layers  of  opposite 
signs.  The  negative  vorticity  layer  is  the  stronger  of  the 
two  and  dominates  the  dynamics.  As  a  result  the 
destabilization  of  the  flow  leads  to  the  formation  of 
downstream  vortical  structures  which  exhibit  substantial 
similarities  to  those  of  the  errorfunction  inlet  profile  flow. 
The  positive  vorticity  is  wrapped  around  the  negative 
vorticity  layer  and  is  finally  entrained  inside  the  vortical 


structures.  Remnants  of  it  survive  far  downstream  in 
small  regions  of  counter-rotating  flow  which  are  very 
beneficial  to  the  local  mixing  but  whose  effect  on  the 
large  scale  dynamics  is  secondary.  The  downstream 
similarity  of  the  errorfunction  and  wake-modified 
flowfields  in  the  uniform  density  case  is  also  evident  in 
the  analysis  of  the  mean  flow.  Mean  velocity  profiles  of 
the  wake-modified  flow,  shown  in  Fig.8  (top)  indicate  that 
within  1.6  channel  heights  downstream  the  wake  deficit 
has  been  annihilated  and  the  velocity  profile  has  been 
altered  to  an  errorfunction  profile,  the  characteristic 
downstream  profile  of  the  errorfunction  flow.  Thus  the 
‘wake  region’,  i.e.  the  region  where  the  mean  flow 
experiences  large  scale  wake  features  is  rather  small.  The 
results  of  Fig.8  (top)  also  provide  an  explanation  for  the 
nonmonotonic  behavior  of  the  shear  region  shown  in 
Fig.7a.  Evidently  the  drop  around  1 .2  <  x  <  1 .6  is  related 
to  the  elimination  of  the  wake  region  and  the 
transformation  of  the  mean  velocity  profile  to  an 
errorfunction. 

The  reacting  field 

The  initial  location  of  the  reaction  interface  with 
respect  to  the  vorticity  layers  is  crucial  to  the  evolution  of 
the  reacting  field.  Tbis  fundamentally  distinguishes  the 
errorfunction  and  wake-modified  flows.  In  the  former, 
the  reaction  interface  is  initially  located  at  the  center  of 
the  single  sign  vorticity  layer,  that  is,  it  coincides  with  the 
surface  of  maximum  (magnitude-wise)  vorticity.  As  this 
layer  gets  destabilized  it  stretches  the  interface  and 
augments  the  burning.  This  behavior  continues  further 
downstream  where  the  instability  of  the  layer  evolves  into 
its  non-linear  stages  and  vortical  structures  form.  In 
contrast,  when  the  wake-modified  inlet  profile  is  used,  the 
reaction  interface  is  located  between  the  two  opposite  sign 
vorticity  layers.  As  the  flow  develops  essenti^ly  through 
the  destabilization  of  the  negative  vorticity  layer,  the 
interface  is  wrapped  around  this  layer  rather  than  being 
inside  it.  This  arrangement  implies  that  the  shear  and 
mixing  regions  of  the  flow  do  not  coincide  initially.  This 
can  clearly  be  seen  in  the  first  two  vortical  structures  of 
the  Go  =  0  wake-modified  flow  in  Figs.  5  and  6  (parts(c), 
top)  where  large  regions  of  shear  do  not  mix  reactants  but 
consist  primarily  of  the  reactant  carried  by  the  fast  stream. 
It  is  only  further  downstream,  where  neighboring  vortical 
structures  pair,  that  the  reaction  interface  is  moved 
towards  the  center  of  the  shear  region  and  its  evolution 
resembles  that  of  the  errorfunction  inlet  profile  flow. 

The  impact  of  these  differences  between  the 
errorfunction  aiKl  the  wake  modified  flows  on  the  burning 
can  clearly  be  assessed  through  the  mean  shear  region  and 
product  integral  plots  of  Fig.7.  For  the  errorfunction  flow, 
the  rise  in  the  amount  of  product  formed  is  closely  related 
to  the  destabilization  of  the  flow  and  the  growth  of  the 
shear  region.  For  the  wake-modified  flow,  on  the  other 
hand,  the  shear  region  and  the  product  integral  do  not 
correlate  initially.  Close  to  the  inlet,  growth  of  the  shear 
region  does  not  translate  to  more  product  Rather  the  rate 
of  product  formation  (judged  by  the  slope  of  the  product 
integral  curve)  is  substantially  increased  around  the  region 
where  the  vorticity  layer  starts  wrapping  the  reaction 
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Velocity  profile,  Wake,  Qo=0,  x=0.8 


Velocity  Profile,  Wake,  Qo=0,  x=1.6 


( 


Velocity  Profile,  Wake,  Qo=6,  x=0.8 


Velocity  Profile,  Wake,  Qo=6,  x=1.6 


Figure  8.  Mean  velocity  profiles  of  the  wake-modified  flow  at  two  downstream  locations,  x  =  0.8 
(left)  and  X  =  1 .6  (right),  for  the  Go  =  0  (top)  and  Go  =  6  (bottom)  cases. 


interface  around  it  (note  approximate  location  in  Figs.5 
and  6.).  Furthermore,  it  continues  at  an  approximately 
constant  rate  fiirther  downstream  as  the  pairing  maintains 
the  mixing  and  burning  at  high  levels.  In  contrast,  the 
shear  region  experiences  a  drop  at  the  early  stages  of  this 
region  of  high  burning.  As  was  shown  earlier  this  drop  is 
related  to  the  elimination  of  the  wake  region  but  here  it 
may  also  serve  as  an  additional  example  of  the  finite 

separation  of  the  shear  and  mixing  regions. 

The  relationship  between  the  shear  and  mixing 
regions  for  the  wake-modified  and  the  errorfunction  flow 
is  further  investigated  in  Fig.9  (top)  and  Fig.  10  (top), 
respectively,  which  contrast  mean  product  mass-fraction 
and  turbulent  kinetic  energy  profiles  at  two  different 
downstream  locations  for  the  wake-modified  and 
errorfunction  flows,  respectively.  The  turbulent  kinetic 

energy  is  defined  as  KE-Q.i.u'^  +  v'^).  The  figures 
clearly  indicate  that  for  the  wake-modified  flow,  the 
cross-stream  locations  of  the  maxima  of  the  kinetic  energy 
and  mass-fraction  profiles  are  not  the  same.  This 


misalignment  is  reduced  downstream  where  the  negative 
vorticity  layer  mixes  the  reactants.  The  profiles 
corresponding  to  the  errorfunction  case,  shown  in  Fig.  10, 
on  the  other  hand,  do  not  experience  any  such 
misalignment  in  the  profiles.  This  implies  that  in  the 
wake-modified  case  turbulence  does  not  necessarily  lead 
to  scalar  mixing  which  is  in  agreement  with  earlier 
statements.  As  will  shortly  be  seen  this  difference 
between  the  errorfunction  and  wake-modified  flows  has 
important  ramifications  when  combustion  exothermicity 
is  present 

V.2.2.  The  exothermic,  variable  density  flow 

The  differences  between  the  errorfunction  and 
the  wake-modified  flows  are  even  more  substantial  in  the 
reacting  case.  This  is  mainly  a  consequence  of  the 
difference  in  the  location  of  the  reaction  interface  with 
respect  to  the  vorticity  field  and  the  resulting 
displacement  of  the  shear  and  mixing  regions  noted 
above.  That  the  over  all  effect  of  exotheimicity  on  the 


L 


10 


k  » 


Figure  9.  Mean  product  mass-fraction  and  turbulent  kinetic  energy  profiles  of  the  wake-modified  flow 
at  two  downstream  locations,  x  =  0.8  (left)  and  x  =  2.4  (right),  for  the  Go  =  0  (top)  and 
Go  =  6  (bottom)  cases. 
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mean  product  formation  is  similar  (i.e.  reduced  in  both 
cases)  is  in  fact  a  result  of  di^erent  mechanisms. 

The  erroifiinctionflow 

In  the  errorfiinction  case  the  deposition  of  energy 
and,  hence,  the  creation  of  volume,  takes  place  at  the 
center  of  the  negative  vorticity  layer  where  the  vorticity 
magnitude  is  maximum.  As  a  result  this  also  becomes  the 
region  of  minimum  density.  This  coincidence  of  the 
vorticity,  density-gradient  and  volumetric  expansion  fields 
has  important  implications  on  the  dynamics  of  the  flow. 
The  reader  is  referred  to  Refs.[4,  6]  where  detailed 
discussions  of  the  vorticity  dynamics  of  a  similar  flow  are 
presented.  Essentially  the  vorticity  field  is  modified  by 
two  major  mechanisms;  volumetric  expansion  and 
baroclinic  generation  (see  Eq.(3)).  The  first  expands  the 
area  and  proportionately  decreases  the  vorticity  to 
preserve  the  circulation.  Consequently,  this  mechanism 


does  not  alter  the  vorticity  sign.  By  thickening  the 
vorticity  layer,  however,  it  delays  the  destabilization  of 
the  flow.  Furthermore,  by  weakening  the  vorticity  it 
results  in  downstream  eddies  with  smaller  rotation  rates 
than  those  of  the  uniform  density  flow.  Both  the  delay  in 
the  onset  of  the  instability  and  the  weakening  of  the 
vorticity  are  major  reasons  for  the  reduction  of  the  shear 
region  and  the  related  (see  earlier  comments)  reduction  in 
the  product  integral  seen  in  Fig.7. 

The  baroclinic  generation  term  is  strongly 
directional  and  modifies  the  vorticity  of  the  fluid  in 
regions  where  the  pressure-gradient  (or  equivalently  the 
fluid  acceleration)  and  the  density-gradient  vectors  are 
misalighed.  The  contribution  of  this  term  tends  to 
improve  the  mixing  within  individual  vortical  structures 
but,  at  the  same  time,  to  inhibit  their  large  scale 
interactions.  This  is  because  this  term  is  the  origin  of  the 
positive  vorticity  surrounding  the  errorfiinction  flow 
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Errorfunction,  Qo=0,  x=1.0 


Errorfunction,  Qo=0,  x=3 


Errorfunction,  Qo=6,  x=l. 


Errorfunction,  Qo=6,  x=3 


YporKE 


Figure  10.  Mean  product  mass-fraction  and  turbulent  kinetic  energy  profiles  of  the  errorfunction  flow 
at  two  downstream  locations,  x  =  1 .0  (left)  and  x  =  3.0  (right),  for  the  Qo  =  0  (top)  and 
j2i>  =  6  (bottom)  cases. 


eddies  seen  in  Fig.5a(bottom)  which  resists  the  clockwise 
pairing  interaction  of  neighboring  structures  imposed  by 
the  free-stream  velocity  ratio.  The  inhibition  of  pairing  is 
another  reason  fix'  the  reduction  in  the  mixing  and  burning 
in  the  errorfunction  flow. 

The  wake-modified  flow 

In  the  wake-modified  case  the  dynamics  are 
quite  different.  Due  to  the  location  of  the  reaction 
interface  the  deposition  of  energy  takes  place  between  the 
two  opposite  sign  vwticity  layers,  i.e.  it  does  not  coincide 
with  the  region  where  vorticity  is  quantitatively  most 
significant.  There  are  two  major  consequences  of  this:  (i) 
The  interaction  of  two  vorticity  layers  is  delayed  as 
conqjared  to  that  of  the  non-exothermic  wake-modified 
flow,  and  (ii)  the  direct  effects  of  volumetric  expansion 
and  baroclinic  generation  on  the  vorticity  dynamics  are 
substantially  diminished  as  compared  to  those  in  the 


exothermic  errorfunction  flow.  These  points  are 
elaborated  upon  in  what  follows. 

The  creation  of  volume  between  the  two  vorticity 
layers  pushes  them  away  from  one  another  and  from  the 
reaction  interface.  As  a  result,  their  interaction  is  delayed. 
This  delay  must  also  be,  to  some  extend,  due  to  the 
weakening  of  the  vorticity  due  to  the  expansion.  Both  of 
these  effects  are  evident  in  Fig.Sc.  For  the  exothermic 
case  the  interaction  of  the  two  layers,  i.e.  the  entrainment 
of  the  positive  vorticity  layer  by  the  negative  vorticity 
layer,  is  substantially  delayed.  The  motion  of  the  layers 
away  from  one  another  close  to  the  inlet  and  the  decay  in 
the  magnitude  of  the  vorticity  there,  are  also  clear.  The 
delayed  interaction  of  the  two  vorticity  layers  essentially 
implies  that  the  wake  region  survives  forther  downstream. 
This  can  be  verified  by  the  mean  velocity  profiles  of  Fig.8 
and  by  the  corresponding  size  of  the  mean  shear  region  of 
Fig.7a.  The  former  indicate  that  the  wake  component  of 
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the  mean  velocity  profile  survives  further  downstream  in 
the  exothermic  case.  The  latter  further  clarifies  that  for 
the  exothermic  case  the  wake  region  does  not  totally  get 
eliminated  within  the  computational  domain. 

Earlier  comments  on  the  importance  of  the  early 
interaction  of  the  two  vorticity  layers  to  combustion 
indicate  the  significance  of  its  delay  in  the  exothermic 
case.  A  reduction  of  mixing  and  burning  is  to  be 
expected,  and  this  was  seen  in  the  product  integral  plots  of 
Fig.7b.  The  instantaneous  flow  visualizations  help 
provide  a  detailed  understanding  of  how  this  comes  about. 
The  lack  of  interaction  between  the  two  vorticity  layers 
leads  to  their  rather  independent  initial  development.  This 
is  particularly  true  of  the  stronger,  negative  vorticity  layer, 
which  rolls  very  much  like  its  errorfunction  counterpart. 
As  in  the  non-exothermic  wake-modified  flow,  however, 
its  effect  on  the  reacting  field  is  fundamentally  different 
from  the  errorfunction  flow  since  the  reaction  interface  is 
no  longer  located  inside  the  negative  vorticity  layer.  In 
contrast  to  the  uniform  density  wake-modifi^  flow,  its 
ability  to  convolute  the  reaction  interface  is  diminished  as 
the  two  are  spaced  further  away  from  one  another  due  to 
the  intervening  volumetric  expansion  filed.  Thus,  vortical 
structures  grow  and  pair  away  from  the  flame  which 
essentially  remains  flat  (see  Figs.  5  and  6).  In  other 
words,  the  shear  mixing  is  initially  irrelevant  to 
combustion.  It  is  only  downstream,  where  the  vortical 
structures  become  strong  enough  that  they  start  wrapping 
the  flame  around  them.  Evidence  of  this  can  be  seen  in 
the  growth  rate  of  the  product  integral  which  experiences 
a  sharp  increase  in  its  growth  rate  when  convolution  of  the 
flame  by  the  vorticity  field  finally  takes  place. 

The  mean  product  mass-fraction  and  turbulent 
kinetic  energy  profiles  of  Figs.9  and  10.  further  clarify 
these  points.  They  indicate  that  in  the  early  sections  of 
the  domain  the  displacement  of  the  shear  and  mixing 
regions  witnessed  in  the  wake-modified  flow  (Fig.9)  is 
enhanced  in  the  presence  of  exothermicity.  As  in  the  non- 
exothermic  flow  this  displacement  is  reduced  downstream 
where  the  negative  vorticity  layer  is  able  to  effectively 
mix  the  reactants.  In  agreement  with  the  non-exothermic 
case  the  errorfunction  does  not  experience  any  such 
displacement  of  the  shear  and  mixing  regions. 

In  the  wake-modified  case,  volumetric  expansion 
occurs  away  from  the  area  where  the  vorticity  magnitude 
is  large.  This  diminishes  its  effect  on  the  vorticity  field. 
As  a  result,  both  of  the  important  consequences  of  the 
presence  of  volumetric  expansion  on  the  errorfunction 
flow,  i.e.  the  delay  in  the  onset  of  the  instability  and  the 
formation  of  weaker  eddies  are  not  as  important  in  the 
wake-modified  flow.  The  instantaneous  visualizations 
indicate  that  the  destabilization  of  the  flow  takes  place  at 
about  the  same  stream  wise  location  for  all  the  wake- 
modified  flow  simulations.  The  effect  of  exothermicity  is 
restricted  to  a  somewhat  lower  growth  rate  of  the 
destabilizing  perturbation.  As  was  already  pointed  out, 
the  reduction  in  the  burning  in  this  case  is  more  closely 
related  to  the  lack  of  convolution  of  the  reaction  interface 
rather  than  to  a  delay  of  the  roll-up  of  the  vorticity  layer 
itself.  As  long  as  the  vorticity  layer  rolls,  sooner  or  later  it 
will  wrap  the  flame  around  it  and  enhance  the  burning. 


This  implies  that  the  wake-modified  flow  is  much  more 
resistant  to  the  negative  impact  of  combustion 
exothermicity  than  the  errorfunction.  (Note  that  the 
specification  of  thicker  inlet  scalar  profiles  for  the 
errorfunction  than  for  the  wake-modified  inlet  condition  - 
see  Section  IV-  implies  that  the  difference  between  the 
errorfunction  and  wake-modified  flows  on  the  issue  of 
instability  suppression  in  the  results  presented  here  is 
actually  underestimated). 

The  displacement  of  the  initial  vorticity  with 
respect  to  the  variable  density  regions  for  the  wake- 
modified  flow  also  diminishes  the  contribution  of  the 
baroclinic  vorticity  generation  term  as  compared  to  that  of 
the  errorfunction  flow.  In  the  latter,  this  mechanism 
resulted  in  the  appearance  of  positive  vorticity  around  the 
vortical  structures.  The  locus  of  the  positive  vorticity  in 
the  visualizations  of  the  wake-modified  flow  shown  in 
Fig.Sc,  on  the  other  hand,  implies  that  this  vorticity  is  a 
result  of  the  inlet  condition  and  not  of  generation.  The 
negative  vorticity  layer  of  the  wake-modified  flow  which 
rolls  like  a  uniform  density  errorfunction  does  not 
experience  any  regions  of  positive  vorticity. 


VL  CONCLUSIONS 

Results  of  two-dimensional  numerical 
simulations  of  spatially-developing  shear  layers  reacting 
at  an  infinite  rate  indicate  that  the  presence  of  a  wake 
deficit  in  the  inlet  velocity  profile  modifies  the 
downstream  dynamics  in  such  a  way  that  mixing  and 
burning  are  increased.  Combustion  exothermicity  reduces 
the  burning  in  this  flow.  While  a  similar  trend  is 
experienced  in  the  flow  arising  from  a  monotonically 
varying  inlet  profile,  in  this  case  an  errorfunction,  the 
mechanisms  by  which  this  is  achieved  in  the  two  flows 
are  different. 

Linear  stability  analysis  of  the  uniform  density 
flow  indicates  that  the  flow  arising  from  the  wake- 
modified  profile  is  much  more  unstable  than  that  from  the 
errorfunction.  This  result  is  confirmed  by  the  numerical 
simulations,  which  also  indicate  that  the  stronger 
instability  of  the  wake-modified  flow  translates  to  better 
mixing  and  increased  burning. 

Beyond  stability  characteristics,  the  most 
important  differences  in  the  behavior  of  the  errorfunction 
and  wake-modified  flows  arise  due  to  the  relative  location 
of  the  reaction  interface  with  respect  to  the  vorticity  field. 
In  the  errorfunction  flow  the  reaction  interface  is  located 
at  the  center  of  a  single,  negative  vorticity  layer,  and  gets 
highly  convoluted  by  it.  For  the  wake-modified  flow  the 
reaction  interface  is  located  between  two  vorticity  layers 
of  opposite  signs.  The  negative  layer,  which  dominates, 
convolutes  the  interface  by  wrapping  it  around  itself. 

The  difference  in  the  location  of  the  reaction 
interface  with  respect  to  the  negative  vorticity  layer 
between  the  errorfunction  and  wake-modified  flows  has 
the  following  important  consequences:  (i)  The  shear  and 
mixing  regions  do  not  always  coincide  in  the  wake- 
modified  flow  as  they  do  in  the  errorfunction.  (ii)  In  the 
wake-modified  flow  exothermicity  decreases  the  product 
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formation  through  a  delay  in  the  convolution  of  the 
reaction  interface  by  the  negative  vorticity  layer  as 
volumetric  expansion  pushes  the  two  away  from  one 
another.  In  the  errorfunction  flow  the  reduction  comes 
about  by  the  weakening  of  the  vorticity  field  by 
volumetric  expansion,  (iii)  The  wake-modified  case  is 
much  more  resistant  to  the  exothermicity  related 
instability  suppression  than  the  errorfunction  case.  This  is 
because  volumetric  expansion  takes  place  away  from  the 
regions  where  the  vorticity  is  most  significant,  (iv)  The 
interaction  of  the  two  vorticity  layers  of  the  wake- 
modified  case  is  delayed  in  the  exothermic  case.  As  a 
result  the  size  of  the  region  is  enlarged  in  this  case. 
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ABSTRACT 

The  objectives  of  this  work  are  to  develop  a 
combustion  submodef  for  high  Damkohler  number 
turbulent  combustion  simulation,  and  to 
investigate  the  transient  dynamics  of  strained  flame. 
A  series  of  mathematical  transformations  are  used, 
under  the  assumption  that  the  applied  strain  rate  is 
constant  throughout  the  flame  structure,  to  translate 
the  system  of  equations  governing  flame 
propagation  in  an  unsteady  stagnation  point  flow 
into  a  set  of  reaction-diffusion  equations. 
Combustion  is  described  by  a  multi-step  chemical 
kinetics  mechanism.  The  flame  structure  computed 
using  this  model  is  compared  with  the  "exact" 
numerical  solutions  obtained  for  unstrained  and 
highly  strained  stoichiometric,  atmospheric 
pressure,  methane-air  flames.  In  the  case  of  an 
unstrained  flame,  our  model  predictions  agree  very 
well  with  the  exact  solution.  In  the  case  of  a  highly 
strained  flame,  the  temperature  and  major  species 
profiles  are  well  predicted,  while  the  minor  species 
( radicals  )  on  the  products  side  are  overpredicted. 
The  transient  dynamics  of  lean  and  stoichiometric 
flames  are  investigated  by  exposing  the  flame  to  a 
sudden  change  in  strain.  For  a  given  value  of  strain, 
we  find  that  the  stoichiometric  flame,  which  has 
the  highest  heat  release  rate,  has  the  smallest 
response  time.  As  the  mixture  becomes  leaner,  the 
heat  release  rate  drops  causing  the  flame  response 
time  to  increase.  For  a  given  equivalence  ratio,  the 
response  time  first  increases  with  strain  reaching  a 
maximum  at  approximately  700  1/s,  then 
monotonically  decreases.  For  strains  lower  than 
700  1/s,  the  flame  is  convected  by  the  flow  while 
its  structure  remains  almost  unchanged  and  the 
flame  dynamics  are  governed  by  the  flame  thermal- 
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diffusion  time  scale  which  increases  with 
equivalence  ratio.  For  strains  higher  than  1000  1/s, 
when  the  influence  of  the  flow  on  the  flame 
structure  is  dominant,  the  flow  time  scale 
determines  the  flame  dynamics  and  makes  it 
independent  of  the  equiv^ence  ratio.  The  results 
show  that  for  turbulent  combustion  simulation,  the 
assumption  of  instantaneous  flame  response  to 
strain  is  accurate  only  for  near-stoichiometric 
flames . 


1.  INTRODUCTION 

Numerical  modeling  and  simulation  of 
turbulent  combustion  in  premixed  gases,  at 
conditions  when  the  flame  thickness  is  smaller  than 
the  smallest  turbulent  eddies,  have  adopted  the 
flamelet  and  thin  flame  approaches,  respectively. 
In  both  approaches,  turbulence,  while  weakly 
changing  the  internal  structure  of  the  flame  and 
thus  its  burning  velocity,  extends  its  overall  surface 
area  and  thus  substantially  increases  the  overall 
burning  rate  per  unit  volume.  This  is  reasonable 
when  the  turbulence  intensity,  as  measured  by, 
e.g.,  the  square  root  of  the  tot^  kinetic  energy  of 
the  fluctuating  field,  is  weak  relative  to  the  burning 
velocity.  As  the  intensity  increases,  local  flow 
gradients  become  comparable  with  the  gradients 
across  the  flame  and  turbulence  induced  strain  rates 
alter  the  flame  structures  in  non  trivial  ways.  These 
effects  are  the  subject  of  this  paper. 

In  our  previous  paper  on  this  subject  [1], 
we  investigated  the  effect  of  the  Lewis  number  and 
the  heat  release  rate  on  the  flame  response  to  a 
sudden  change  in  the  strain  rate  and  to  an 
oscillating  strain  using  a  “uniform”  strain  model. 
We  have  shown  that  both  the  heat  release  and  the 
deviation  of  the  Lewis  number  from  unity  can  have 
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profound  effects  on  the  flame  dynamics.  A  simple 
one-step  chemical  mechanism  and  constant 
thermodynamic  properties  were  used  in  these 
studies.  While  providing  insight  into  the 
mechanisms  governing  the  flame  dynamics,  that 
model  is  not  capable  of  predicting  the  dynamics  of 
real  hydrocarbon-air  flames  due  to  the  simplicity  of 
the  chemical  kinetics  and  transport  mechanisms 
utilized. 


The  objectives  of  this  work  are  to 
investigate  the  performance  of  our  model  when  used 
in  situations  where  multi-step  kinetics  are 
important,  and  to  determine  the  flame  response 
time  when  unsteady  perturbations  are  imposed.  We 
show  that  our  model,  which  makes  the 
simplification  that  the  strain  rate  is  constant 
throughout  the  flame  structure,  predicts  accurately 
the  flame  structure  in  terms  of  the  major  species, 
over  a  wide  range  of  strain  rates.  This 
simplification  allows  one  to  transform  the 
governing  equations  into  reaction-diffusion 
equations  and  thus  reduces  the  computational  effort 
required  to  obtain  the  solution,  especially  in  the 
unsteady  case.  We  also  find  that  the  response  time 
of  the  flame  is  a  strong  function  of  the  heat  release 
rate  and  the  applied  strain. 

This  paper  represents  a  logical  application 
of  the  approach  developed  in  [1]  to  flames 
described  by  the  multi-step  chemical  kinetics.  The 
problem  of  flame  propagation  in  mixtures  with  a 
multi-step  kinetics  mechanism  has  been  studied 
before  using  different  models  (  see,  for  example, 
[2],  [5],  [7],  [8]-[12]  among  others).  Most  of  these 
studies  assume  steady  flames  and,  Aus,  may  not  be 
applicable  as  sub-models  in  a  turbulent  combustion 
simulation  in  situations  when  the  dynamics  of 
flame-flow  interactions  are  important.  As  we 
showed  in  our  previous  studies  [1],  [13],  flame 
response  to  unsteady  perturbation  can  have  a 
significant  impact  on  the  burning  rate,  extinction 
and  combustion  instability. 

The  paper  is  organized  as  follows.  In 
Section  2,  the  unsteady  problem  is  formulated  and 
the  solution  procedure  is  described.  Section  3  deals 
with  the  numerical  aspects  of  the  problem.  In 
Section  4,  we  discuss  the  results  of  the 
calculations,  compare  the  model  with  two  exact 
numerical  solutions  for  the  cases  of  unstrained  and 
strained  flames  and  show  the  effect  of  a  step-wise 
variation  in  the  strain  on  the  response  time  at 
different  equivalence  ratios.  Section  5  contains  the 
conclusions. 


2.1  FORMULATION 

The  flow  configuration  is  depicted  in 
figure  1.  The  flame  is  located  in  the  stagnation 
point  flow  produced  by  two  counter-flowing  planar 
jets  of  reactants  and  pixxlucts.  The  x-  and  y-axes  are 
directed  parallel  and  perpendicular  to  the  flame, 
respectively.  The  flame  is,  in  general,  unsteady, 
i.e.  its  location  and  burning  velocity  change  with 
time  due  to  changes  in  the  flow  field,  species 
concentration  ot  the  temperature. 

Applying  the  boundary  layer 
approximation,  in  which  diffusion  along  the 
boundary  layer,  taken  here  to  be  the  x-axis,  is 
negligibly  small,  we  write  the  equations  governing 
the  flame  structure  as  follows: 


dt  dx  dy 


(2.1) 


+  pv  ^=.±{pYkVk)  + 

ax  ay  ay 

6kWk,fork  =  0,K  (2.2) 


dt  dx  dy 


dy\  dyl 


-  X  pYkVkCp. k^-Z(^kHk,  (2.3) 
dy 


p=plu-T,  (2.4) 

Wmix 

p^  +  p«fi  +  pv^  =  -f,  (2.5) 

dt  dx  dy  dx 


where  p  is  the  mixture  density;  u  and  v  are  the 
velocities  of  fluid  in  the  x  and  y  directions, 
respectively;  7*  is  the  mass  fraction;  D*  is  the 
mixture  averaged  diffusion  coefficient;  m  is  the 
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production  rate;  IV*  is  the  molecular  weight;  Cp  and 
Cp,  k  are  the  constant  pressure  specific  heats  of  the 
mixture  and  individual  species,  respectively;  (K+1) 
is  the  total  number  of  chemical  species;  T  and  p 
are  the  thermodynamic  temperature  and  pressure, 
respectively;  A  is  the  thermal  conductivity  of  the 
mixture;  /?«  is  the  universal  gas  constant;  is 
the  molecular  weight  of  the  mixture;  Hk  is  the 
enthalpy  in  molar  units;  and  the  diffusion  velocity 
Vk  is  defined  as  follows: 


Xk  dy  Xk  T  dy 


for  k  =  0,K 


P  =P 


Ru 


T  . 


(2.9) 


while  the  x-momentum  equation  is  redundant  along 
the  stagnation  streamline.  The  latter  simplification 
makes  our  model  different  from  the  models 
described  in  Refs.  [2],  [5],  [7],  [8]  -  [12],  where  the 
momentum  equation  is  solv^.  The  implication  of 
not  satisfying  the  momentum  equation  will  become 
clear  when  we  compare  the  model  results  with  the 
exact  solution. 


where  Vcor  is  the  correction  velocity  introduced  in 

K 

order  to  satisfy  the  identity  Y  YkVk  -0  ;  kr.k  is 

k=0 

the  thermal  diffusion  coefficient.  In  the  following, 
we  neglect  the  thermal  diffusion  velocity  (the  last 
term  of  the  above  equation)  and  the  correction 
velocity.  As  shown  in  [2],  this  simplification  does 
not  significantly  affect  the  results.  In  the  above 
equations,  we  neglected  the  diffusion  terms  in  the  x 
-momentum  equation,  and  the  heat  and  mass 
transport  equations,  and  the  dissipation  term  in  the 
energy  equation. 

Along  the  stagnation  streamline  x  =  0, 

and  due  to  symmetry  :  u  =  0  ,  —  =  0,  leading  to 

dx 

the  following  simplifications  of  equations  (2.1-5): 


The  boundary  conditions  for  equations  (2.6-8)  are: 

;y  =  .oo  ;  Yk  =  Yk.b.  T  =  Tb  ,  0  (2.10) 

dy 

y  =  +oo  :  Yk  =  Yk.u,  T  =  Tu  (2.11) 

where  the  subscripts  u  and  b  denote  values  on  the 
unburnt  and  the  burnt  side,  respectively.  The 
boundary  conditions  for  the  velocity  field  will  be 
discussed  later.  Although  along  the  stagnation 
streamline  u  =  0,  the  ^  !  dx  term  of  equation 
(2.6)  does  not  vanish.  Since  the  boundary 
conditions  for  T*  and  T  are  independent  of  x,  it  can 
be  demonstrated  [1]  that  the  flame  structure  near  the 
stagnation  streamline  is  similar  to  that  described  by 
equations  (2.6-9). 


^+p^+^(M.  =  o  (2.6) 

dt  dx  dy 

p  ^  +  p  V  d(W^ixYk)\  ^ 

dt  dy  dy  dy  I 

(6kWk,k  =  0,K  (2.7) 


One  more  simplification  of  the  problem  is 
made.  The  energy  equation  (2.8)  contains  a  term  in 
the  form 

i  pYkVkCp.k^. 
k=o  dy 

Smooke  [2]  showed  that  this  term  can  be  safely 
neglected  virtually  for  all  important  configurations 
of  strained  and  unstrained  flmes.  Eliminating  this 
term,  equation  (2.1 1)  becomes: 


PCp 


dt 


K- 

d 

(,  dr\ 

pcp(  — 

.V^)  = 

d_ 

(^-1 

dy 

dy 

H)- 

dt 

dy 

dy 

1  dyl 

K  dr  ^  ■ 

2^  pYkVkCp,  k  0)k  Hk  (2.8) 


k=0 


k^O 


(2.12) 
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The  solution  technique  applied  here  to  solve  the 
system  of  equations  (2.6-12)  with  boundary 
conditions  (2.10-11)  is  identical  to  that  used  in  our 
previous  work  [1],  Hence,  we  will  describe  it  only 
briefly.  First,  we  decompose  the  velocity  field  as 
follows: 


V  =  V  +  Vr 


(2.13) 


where  v  is  the  total  velocity,  is  the  combustion 

generated  velocity,  and  V  is  the  externally  imposed 
velocity  which  satisfies  the  following 
incompressible  continuity  equation: 


(2.14) 


Substituting  equation  (2.13)  into  continuity 
equation  (2.6),  using  equation  (2.14)  and  assuming 
that  on  the  stagnation  streamline  duc/  dx  =  0,we 
obtain: 


+v  ^  =  0  (2.15) 

dt  dy  dy 

Substitution  the  velocity  decomposition  in 
equation  (2.13)  into  equations  (2.7)  and  (2.12) 
yields: 


ct 


y  =  exp( 


e(t')dt')  y  ^B(t)y,  t  =  t  (2.19) 


0 


eliminates  the  V  d!  dy  operators  from  equations 
(2.15-17).  Following  the  application  of  the 
transformation  in  equation  (2.19),  the  governing 
equations  take  the  form: 


=  (2.20) 

A  ^ 


P 


dt 

b2±1  pDk 

^  Wmix 


dYk_ 

d(W,nixYk)\ 

dy  I 


+ 


m  Wk,  k  =0.  K 

(2.21) 


+  pCp  VcB 


B^— 

dy 


X  ^klik 


k=0 


(2.22) 


P  ^  +  P  (V+Vc) 
dt  dy 

d  /  pDk  djWmixYk) 

dy  Wfft/jc  dy 


+  (bkWk,k=0,K 
(2.16) 


+  pCp(V+Vc)  ^  = 

dy 


M)-  i 


cOkHk 


(2.17) 


Second,  if  the  externally  imposed  inviscid 

— * 

irrotational  velocity  field  V  is  specified  as 

V=(e(t)x,-e(t)y)  (2.18) 

where  e  is  the  strain  rate,  then  the  following 
U-ansformation  [3] 


The  combustion  generated  velocity  is 

eliminated  by  applying  another  (Howarth) 
transformation  to  equations  (2.20-22): 


(2.23) 


where  the  oveibar  denotes  the  density  normalized  by 
its  value  on  the  reactants  side,  pu.  To  apply  this 
transformation,  we  divide  the  continuity  equation 
by  pu  and  integrate  it  to  obtain: 


^  dy+  B  p(y)  vc(y)  -  B  p(0)  Vc(0)  =  0 

J  dt 

J  y 


(2.24) 
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At  the  stagnation  point  y  =  0  the  total  velocity 
must  be  zero.  This  requirement  imposes  the 
following  boundary  condition  Vc(0)  =  0  and 
equation  (2.24)  is  reduced  to: 


^dy  =-Bp(y)vc(y) 
dt 


(2.25) 


Mapping  equations  (2.21-22)  onto  the  (^,t) 
domain,  and  using  the  chain  rule  and  equation 
(2.25),  we  get: 


O  ^  -j2-  ^  pDk  pd(WmixYk) 
dr  d^Wmix 

+  d>kWk  fork  =  0.K  (2.26) 

pc^=B2p^Up^\-  t<OkHk  (2.27) 
dr  dq  \  dql  k=o 


Finally,  equations  (2.26-27)  are  non- 
dimensionalized  using  the  length  sc^e  based  on  the 

thermal  diffusion  length  y^cl  =  Va^  tscU  where 
Uu  is  the  thermal  diffusivity  coefficient  of  the 
mixture  on  the  reactants  side,  tsci  =  ««  ISu  ,  and  5« 
is  the  unstrained  flame  burning  velocity. 
Introducing  the  non-dimensional  temperature, 
9  =(T  -Tu)/  (Tb  -  Tu),  we  get: 


p  ^  =B^p  jL^ij^pi(Ern^] 

dr  CCu  \Wmix  / 


+  tscimWk  fork  =  0,K  (2.28) 


=  B^  p 

a« 


hcl  V 
(  Tb’Tu)  k^o 


mHk, 


(2.29) 


where  ^  =  ^  Jysci ,  and  r  =  rl  tscu  and  the 
boundary  conditions  are  given  by  equations  (2.10- 
11). 


Next,  we  discuss  the  evaluation  of  the 
transport  coefficients  A  and  D^,  k  =  1,K.  We  used 
the  SANDIA  Transport  [4]  and  Premixed  Flame  [5] 
codes  to  study  the  variation  of  these  parameters 
with  temperature.  In  accordance  with  Smooke  [2], 
we  found  that  the  Lewis  numbers 
^e,k=  ^  (P  Cp  ^k )  of  all  species  in  the 
stoichiometric,  atmospheric  pressure,  methane/  air 
flame  are  almost  constant  with  temperature. 
Therefore,  in  our  calculations  we  used,  the 
following  values  for  the  Lewis  numbers: 

Table  1.  Constant  Lewis  Numbers  of 
Different  Species. 


Species  Lewis  Number 

CN4 . 0.97 

02 . 1.11 

ff20 . 0.83 

CO2 . 1.39 

H . 0.18 

O . 0.70 

OH . 0.73 

HO2 . 1.10 

H2 . 0.30 

CO . 1.10 

H2O2 . 1.12 

HCO . 1.27 

CH2O . 1.28 

C//3 . 1.00 

CH2O . 1.30 

N2 . 1.00 


For  an  accurate  simulation  it  is  extremely 
important  to  have  the  correct  dependence  of  the 
thermal  conductivity  coefficient  on  temperature.  At 
the  same  time,  its  direct  calculations  using  the 
collisions  integrals  could  be  computationally 
prohibitive.  Using  the  SANDIA  Transport  code,  we 
found  that  the  interpolation  for  the  thermal 
conductivity  coefficient  suggested  by  Smooke  [2] , 

X  =  2.5S  10-^ Cp  (Tl  Tu)^  '^  (CGS  units),  produces 
results  which  are  virtually  identical  to  those 
obtained  from  the  code  over  the  entire  range  of 
temperature.  Thus,  this  interpolation  formula  has 
been  adopted  in  our  model. 

Using  the  assumption  of  constant  Lewis 
number,  Lg,  t  =  XI  (p  CpD^)  =  const,  and  the 

Smooke  correlation  X  =  258  W^Cp  (T/ 
we  rewrite  the  diffusion  terms  in  the  governing 
equations  (2.28-29).  For  the  mass  fraction  equation: 
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pDiP=  J«-= 

Cp  k  k  ^u  ^P  ^P>  “ 

-l-.J^SPdLp  au  ,  (2.30) 

^e,  k  ^P 

where  we  used  the  definition  of  the  thermal 
diffusivity  coefficient,  au  =  ^u/  Pu  Cp,  u  •  For 
the  energy  equation: 

Ap=  X  f  =  X^^^SEdL  =  auP  Cp, 

Pu  Au  Pu  Cp,  u 

(2.31) 

Using  Smooke's  correlation  and  the  equation  of 
state,  we  evaluate: 


X  _ 

t' 

0.7 

1-  = 

Cp 

t' 

Xu  <^p 

Ju. 

> 

Xu 

Cp,u 

Ju. 

the  given  mass  fraction  vector  [Yo,  Yu,..,,YkY 
and  temperature  are  evaluated  as  follows: 


=  I 


f  Cp.  k  Yk 


Wk 


^mix  - 


y  Yjl 

^0  Wk  J 


P  =  £jw  (2.35) 

RuT 

where  Cp  is  the  specific  heat  of  mixture  in  mass 
units. 


The  following  mechanism  of  18 
elementary  non-reversible  reactions  and  13  reacting 
species  was  used  to  describe  the  methane-air 
atmospheric  pressure  combustion  [7]: 

Table  2.  The  chemical  kinetics  mechanism  used 
in  the  model. 

SPECIES  : 

CH4  CH3  CH20  HCO  C02  CO 
H2  H  02  O  OH  H02  H20  N2# 


P  _  ^ mix  Yfi  (2.32) 

Pu  W,nix,  u  T 

Dividing  equations  (2.28-29)  by  p  and  p  Cp, 
respectively,  and  using  equations  (2.32),  we  obtain: 


—  =  B‘ 
<5? 


l^e.  k  1 

^scl 

p 

’■  ±Al 

^p  d^\ 

it 

hcl 

j 

y^mix 

au 

Wmix,  u 

Wi 

p  Cp(  Tb-Tu)  *=0 


flit  Hk 


A  ( .Emx..  Yk)\  + 


(2.33) 


(2.34) 


with  boundary  conditions  given  by  equations  (2.10- 
2.11). 


CH4+H=>CH3+H2 

;  2.2E4 

3.0 

8527.7 

CH4+0H=>CH3+H20 

;  1.6E6 

2.1 

2452.2 

CH3+0=>CH20+H 

;  7.0E13 

0.0 

0.0 

CH20+H=>HC0+H2 

:  2.5E13 

0.0 

40003.3 

CH20+0H=>HC0+H20 ;  3.0E13 

1  0.0 

1199.8 

HCO+H=>CO+H2 

;  2.0E14 

0.0 

0.0 

HCO+M=>CO+H+M  ; 

7.14E14 

0.0 

16811.7 

HC0+02=>C0+H02  ; 

3.0E12 

0.0 

0.0 

CO+OH=>C02+H 

;  4.4E6 

1.5 

-740.9 

C02+H=>C0+0H  ; 

1.6E14 

0.0 

26338.4 

H+02=>0H+0  ; 

1.2E17 

1 

O 

’o 

16532.0 

0H+0=>H+02 

:  1.8E13 

0.0 

0.0 

0+H2=>0H+H 

;  1.5E7 

2.0 

7564.5 

0H+H=>0+H2  ; 

6.7E6 

2.0 

5573.6 

0H+H2=>H20+H 

;  1.0E8 

1.6 

3303.0 

H+H20=>0H+H2  ; 

4.6E8 

1.6 

18601.8 

2*0H=>H20+0  ; 

1.5E9 

1.14 

0.0 

H20+0=>2*0H  ; 

1.5E10 

1.14 

17282.5 

H+02+M=>H02+M 

;  2.0E18 

-0.8 

0.0 

H02+M=>H+02+M  ; 

;  2.77E18 

-0.8 

46701.7 

H+0H+M=>H20+M 

:  2.15E22 

-2.0 

0.0 

H20+M=>H+OH+M  ; 

3.87E23-2.0  : 

119397.7 

H02+H=>2*0H  : 

1.5E14 

0.0 

1001.4 

H02+H=>H2+02 

:  2.5E13 

0.0 

690.7 

H02+0H=>H20+02 

;  2.0E13 

0.0 

0.0 

The  standard- state  thermodynamic 
properties  are  evaluated  in  terms  of  the  polynomial 
fits  of  the  specific  heats  at  constant  pressure,  as 
done  in  Chemkin  II  [6].  The  coefficients  of  the  fits 
are  taken  from  the  SANDIA  thermal  data  base  file 
"THERMDAT".  The  properties  of  the  mixture  for 


where  each  line  describing  a  reaction  contains  the 
reaction  equation,  the  pre-exponential  factor,  the 
temperature  exponent,  and  the  activation  energy  of 
the  forward  reaction  constant,  respectively.  All 
numerical  values  are  in  CGS  units  except  for  the 
activation  energy  which  is  given  in  cal  /  mole. 
This  mechanism  is  capable  of  predicting  with 
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reasonable  accuracy  the  burning  velocity  and  the 
profiles  of  temperature  and  species  mass  fractions 
in  lean  to  stoichiometric  flames.  The  mechanism  is 
read  by  an  interpreter  function  which  sets  up  the 
chemical  kinetics  and  thermal  data  bases.  The 
format  of  the  mechanism  description  file  is  very 
close  to  that  used  in  Chemkin  II  [6].  The 
calculation  of  the  species  production  rates  is  also 
very  similar  and  will  not  be  described  here.  For  the 
above  reaction  mechanism,  the  values  of  the  third 
body  efficiencies  are  give  in  [7]  as 


=  h  (XH2,i=L  002,1  =  0.4 

OH2O,  i  =  65  ,  aco,  i  =  0.75  ,  aco2,  i  = 

OCH4 ,  i  =  ^-54,  aN2.  i  =  0.4. 

The  values  of  third-body  efficiencies  for  other 
species  are  equal  to  zero. 


3.  NUMERICAL  SOLUTION 


The  numerical  solution  of  equations  (2.33- 
34)  with  boundary  condition  (2.10-11)  is  obtained 
using  a  Crank-Nicholson  finite-difference  scheme 
using  an  implicit  approximation  of  the  source  term 
to  maintain  stability,  and  the  Thomas  tri-diagonal 
matrix  inversion  dgorithm.  The  details  of  the 
numerical  procedure  were  described  extensively  in 
our  previous  work  [1]  where  a  solution  for  the  one- 
step  chemistry  model  was  obtained.  Here  we 
present  only  the  features  of  the  numerical  procedure 
which  are  different  when  a  multi-step  chemical 
kinetics  mechanism  is  considered. 

First,  in  the  case  of  multi-step  ch^istry, 

the  grid  in  the  computational  domain  (t,  |)  must 
be  nonuniform  because  the  reaction  zone  is  thinner 
than  the  overall  flame  structure  thickness.  This  is 
due  to  the  difference  in  the  time  scales  of  the  fast 
fuel  consumption  and  the  slow  C02  formation 
reactions.  Following  a  rapid  initial  increase  in 
temperature  in  the  reaction  zone,  a  relatively  long 
zone,  of  the  order  of  10  cm  under  normal 
conditions,  is  required  for  the  temperature  to  reach 
its  equilibrium  value.  Every  time  step,  the  grid  is 
adapted  in  such  a  way  that  the  regions  of  strong 
gradients  of  the  dependent  variables 


Yq,  Y],  ....  Yk,  6  are  well  resolved  following  the 
criteria : 


iTi+y  -fi  l<  «  (f)  -  fnin(f)  I  (3.1) 

P  I  max  -  min  {^/^Dl 

where  a  and  jS  are  the  adaptation  tolerance 
parameters  of  the  order  of  0.1;  fi  and  {dfl  are 
the  values  of  a  dependent  variable  and  its  derivative 
at  grid  point  i .  If  criteria  (3.1)  are  not  satisfied,  a 
new  grid  point  is  introduced  half  way  between  the 
points /+7  and  I.  The  values  of  mass  fractions  and 
temperature  at  this  midpoint  are  calculated  as  the 

average  values:  fi+1/2  =  ^(fi+l+fi)  •  This 

interpolation  maintains  the  second  order  of 
accuracy. 


Second,  at  the  onset  of  every  time  step 
the  thermal  thickness  of  the  flame  is  compared  to 
the  distance  from  the  flame  to  the  upstream  and 
downstream  boundaries  of  the  computational 
domain.  If  either  of  these  distances  is  smaller  than 
two  flame  thicknesses,  the  computational  domain 
is  legridded  by  moving  the  boundaries  further  away 
from  the  flame.  The  values  of  the  mass  fractions 
and  the  temperature  at  the  added  mesh  points  are 
taken  to  be  equal  to  the  boundary  values  of  the 
"old"  mesh.  Hence,  the  zero-gradient  boundary 
conditions  are  naturally  satisfied  on  the  new  mesh. 
Each  time  step,  the  physical  coordinate  y,-  is 
recalculated  in  such  a  way  that  the  stagnation  point, 
yi  =  0 ,  remains  fixed. 

Third,  the  time  step  of  the  calculation  is 
taken  to  be  the  minimum  of  the  characteristic  time 
scales  of  species  production: 


hcl,  k  — 


j  V^k,  max3 
3  max 


(3.2) 


where  max  {oflfe}  is  the  maximum  of  the 
production  rate  of  the  k  -th  species,  [Xjt,  maxi  is  its 
molar  concentration  at  this  maximum;  and  the 
computational  time  scale  based  on  a  Courant 
criterion: 

hcl,  Courant  —  twin  {z\  exp  {  -  2  Kq  X  ),  (3.3) 
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i.e. 

=  tnin  {  tgcl,  Of  hcl,  hhcl,  2>- 

•fhcU  Kfhcl,  Courant^  2  j  .  (3.4) 

Equation  (3.4)  shows  that  the  Courant  time  step  is 
an  exponentially  decreasing  function  of  time.  For 
highly  strained  flames,  at  large  Ka,  the  time  step 
decreases  rapidly  with  time,  exponentially  slowing 
down  the  calculation.  A  way  to  avoid  this  problem 
is  to  regrid  the  computational  domain  in  such  a 

way  that  min  (a  ^  is  increased  whenever  the 
time  step  becomes  exceedingly  small. 

We  use  the  profiles  provided  by  the 
SANDIA  Premix  Code  with  the  chemical  kinetics 
model  described  in  Table  2  as  the  initial  profiles. 


4.  RESULTS  AND  DISCUSSION 
Comparison  with  Test  Profiles 

To  establish  the  validity  of  our  model,  we 
tested  the  steady-state  species  mass  fraction  and 
temperature  profiles  obtained  by  our  model  in  the 
case  of  unstrained  and  strained  flame  against  those 
of  Smooke  [2]  and  Rogg  [7],  respectively.  In  the 
case  of  the  unstrained  flame,  the  profiles  were 
plotted  against  the  following  normalized  coordinate: 

y  =  Pu  Su  (y- jf)  cpf  ^  (4.1) 

In  figure  2  we  show  the  temperature  and 
the  reactants  mole  fraction  profiles,  CH4  and  02- 
The  agreement  between  the  model  prediction  and 
the  exact  solution  is  good,  although  the  model 
slightly  overpredicts  the  temperature  in  the  flame 
region,  while  it  underpredicts  it  on  the  reactants 
side.  Consistent  with  these  temperature  profiles, 
our  model  produces  slightly  lower  values  for  Xcha 
behind  the  reaction  zone,  y  =  0  ,  than  these 
obtained  in  Smooke's  model  [2] 

In  figures  3-4  the  major  species  H2O  , 
H2,  CO,  and  CO2,  and  the  radical  H  mole  fraction 
profiles  are  plotted.  The  agreement  between  the 
model  results  and  the  exact  solution  is  very  good. 
The  strongest  disagreement  is  detected  in  the  values 
of  H.  Although  the  part  of  the  profile  on  the 
reactants  side  and  its  maximum  value  are  actually 
captured  by  our  model,  the  post-flame  values  are 


lower  than  that  predicted  in  [2].  This  could  be 
explained  by  the  strong  sensitivity  of  the  H  profile 
to  the  temperature  which  is  overpredicted  by  our 
model  in  the  post-flame  region  (see  figure  2).  Since 
our  model  predicts  accurately  the  major  species 
profiles,  we  expect  the  burning  velocity  to  be  close 
to  the  exact  value.  The  value  provided  by  our 
model,  38  cm/s,  is  very  close  indeed  to  the  exact 
solution  value  of  37.67  cm/s  [2]. 

In  the  case  of  a  highly  strained  flame,  the 
test  profiles  of  Rogg  [7]  and  our  profiles  are  plotted 
in  broken  and  solid  lines  in  figures  5-8.  The 
profiles  of  Rogg  [7]  wCTe  shifted  in  such  a  way  that 
the  positions  of  the  fuel  consumption  layers  in 
both  models  match.  The  temperature  profiles  are 
shown  in  figure  5.  Our  model’s  profiles  and  those 
obtained  using  the  original  equations  exhibit  very 
close  similarity,  especially  in  the  reaction  zone 
region.  In  the  post-flame  zone,  the  temperature  is 
higher  in  Rogg’s  results  [7].  It  should  be  noted  that 
the  equilibrium  values  of  the  temperature  in  both 
models  are  the  same.  However,  the  part  of  the 
domain  shown  in  the  figure  is  not  large  enough  to 
include  the  "hot"  boundary  used  in  the  calculations. 
The  exact  and  model  results  for  CH4 ,  H2O  and  O2 
mass  fractions  are  virtually  identical.  The  CO2  and 
CO  mass  fraction  profiles  are  shown  in  figure  6. 
The  profiles  are  identical  in  the  reaction  zone  but 
differ  slightly  in  the  post-flame  region.  The  profile 
of  CO2  predicted  in  [7]  has  slightly  higher  values 
of  CO2  mass  fraction  here.  The  same  can  be  said 
about  the  CO  mass  fraction  profiles.  Clearly,  the 
model  predicts  well  the  temperature  and  major 
species  profiles  with  the  maximum  discrepancy  in 
the  CO  profile  being  close  to  10  %  . 

A  more  stringent  test  for  the  model  is  how 
accurately  it  predicts  the  minor  species  profiles. 
These  species  are  present  in  very  low 
concentrations,  and  many  of  them  are  extremely 
sensitive  to  the  temperature  and  chemical  kinetics 
mechanism.  One  of  these  species  is  HO2  ,  which  is 
shown  in  figure  7.  The  agreement  between  the 
model  predictions  and  the  exact  solution  is  good, 
although  the  peak  value  is  slightly  overpredicted.  In 
figure  8,  OH,  H,  and  0  radicals  profiles  are 
depicted.  Within  the  reaction  zone,  our  profiles 
match  well  the  exact  solution.  The  post-flame 
values  are,  however,  overpredicted  by  our  model. 

The  above  indicates  that  our  model 
accurately  predicts  the  integral  characteristics  of  the 
flame,  e.g.  the  burning  velocity,  and  the  major 
species  profiles  over  a  wide  range  of  strains. 
Moreover,  at  all  values  of  strain,  the  minor  species 
are  well  predicted  within  the  reaction  zone  but 
deviate  somewhat  on  the  products  side.  The  reason 
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for  the  latter  appears  to  be  the  absence  of  the 
momentum  equation  in  our  model.  Although  the 
momentum  equation  is  redundant  along  the 
stagnation  streamline,  it  provides  means  by  which 
the  condition 

€b  =  ^  Pu  Ipb  » 


Q->  ai-¥a2exp(‘tl  X), 


where  Q  ^ 


G)kHkdy  , 


which  relates  the  strains  on  the  reactants  and 
products  sides,  can  be  imposed.  This  condition 
implies  that  for  the  stoichiometric,  atmospheric 
pressure  flame  the  strain  on  the  products  side  is 
approximately  three  times  larger  than  that  on  the 
reaction  side.  Our  model  does  not  satisfy  this 
boundary  condition,  accounting  for  the  density 
variation  in  the  reaction  zone  by  the  combustion 
velocity  only,  and  hence  it  is  does  not  allow  the 
strain  to  vary  throughout  the  flame  structure.  Thus, 
effectively,  while  the  flame  structure  is  resolved 
correctly,  the  strain  rate  on  the  products  side  in  our 
model  is  lower  than  its  actual  value. 

The  Transient  Dynamics  of  Strained  Flames 

The  response  of  the  flame  to  a  step-wise 
change  in  the  strain  rate  is  now  investigated  using 
results  of  lean  and  stoichiometric  flames.  To  vary 
the  equivalence  ratio,  the  number  of  moles  of 
oxygen  and  methane  are  changed  in  such  a  way  that 
the  total  number  of  the  moles  remains  constant. 
Three  values  of  the  equivalence  ratio  are  considered: 
0  =  0.6,  0.5,  7.0,  To  characterize  these  flames,  the 
steady-state  values  of  the  total  heat  release  rate  are 
plotted  in  figure  9  as  a  function  of  the  strain  rate 
and  equivalence  ratio  .  For  the  unstrained 
stoichiometric  flame  this  value  is  close  to  0.1 
kJ!  cm^-s.  As  the  strain  rate  increases,  the  value 
gradually  decreases  to  0.07  kJi  cm^-s  at  e  =  1800 
1/s,  corresponding  to  the  Karlovitz  number  of  0.18. 
For  leaner  mixtures  with  equivalence  ratios  of  0.8 
and  0.6,  the  total  heat  release  is  approximately  0.08 
and  0.04  kJ!  cm^-s,  respectively.  Again,  the  heat 
release  rate  decreases  monotonically  as  the  strain 
rate  increases  and  the  flame  £q)proaches  extinction. 
Due  to  the  boundary  conditions  (2.10-11),  however, 
extinction  is  only  an  asymptotic  state. 

We  found  that  the  total  heat  release 
determines  the  flame  dynamic  response.  A  typical 
response  of  a  flame  with  multi-step  kinetics  to  a 
sudden  change  in  the  strain  rate,  shown  in  figure 
10,  is  very  similar  to  that  observed  when  using  a 
one-step  Wnetics.  To  quantify  the  flame  response 
time,  the  total  heat  release  as  a  function  of  time  is 
interpolated  using  the  following  exponential 
function 


and  the  coefficients  aj,  a2  and  t  are  calculated 
using  a  least-square  interpolation  procedure. 

In  figure  1 1  the  flame  response  time,  t  , 
is  shown  as  a  function  of  the  equivalence  ratio  and 
the  applied  strain.  For  every  equivalence  ratio 
considered,  the  time  constant  increases  with  strain, 
reaching  a  maximum  around  700  1/s  and  then 
decreases.  This  value  of  strain  corresponds  to 
Ka  =  0.07,  0.103,  0.356  for  0  =7,  0.5,  0.6, 
respectively,  .  The  impact  of  the  strain  on  the 
response  time  can  be  explained  as  follows.  In  our 
previous  paper  [1],  where  a  one-step  kinetics  model 
was  used,  we  show  that  an  increase  in  the  heat 
release  rate  makes  the  flame  responded  faster  to 
strain.  Consistent  with  this  result,  the 
stoichiometric  flame  which  has  the  highest  heat 
release  rate  has  the  smallest  response  time.  The 
maximum  value  of  the  t  (  0,  £  )  curve  can  be 
explained  by  the  competition  between  the  different 
mechanisms  controlling  the  flame  dynamics.  For 
lower  values  of  strain,  the  influence  of  the  flow 
manifests  itself  by  moving  the  entire  flame 
structure  to  a  location  in  the  stagnation  flow  where 
the  mass  conservation  principle  is  satisfied.  Within 
this  range  of  strain  rate,  the  structure  of  the  flame 
does  not  change  and  the  flame  dynamics  is 
controlled  primarily  by  the  thermal-diffusion  time 
scale,  which  is  JO-^.  1.47  5.11  lO"  s. 

for0  =1,0. 8, 0.6,  respectively.  As  figure  11 
indicates,  the  switching  from  the  regime  controlled 
by  the  flame  time  scale  to  the  regime  controlled  by 
the  flow  occurs  around  700  1/s.  For  higher  strains, 
the  influence  of  the  flow  on  the  flame  becomes  so 
overwhelming  that  the  flow  time  scale,  1/  e,  starts 
to  control  the  flame  response.  When  the  strain  rate 
increases,  the  values  of  the  response  times  of 
flames  with  different  equivalence  ratios  become 
closer  to  each  other  and  t  ( <!>,  e )  curves 
converge. 


5.  CONCLUSIONS 

1.  A  model  of  unsteady  strained  flame 
with  multi-step  chemical  kinetics  has  been 
developed  by  applying  a  series  of  mathematical 
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transformations  to  the  system  of  equations 
governing  flame  propagation  in  a  stagnation  flow, 
and  translating  the  system  into  a  set  of  reaction- 
diffusion  equations.  The  flame  structure  calculated 
by  the  model  were  compared  with  the  exact 
solution  for  the  cases  of  unstrained  and  highly 
strained  stoichiometric,  atmospheric  pressure, 
methane-air  flames.  In  the  case  of  unstrained  flame, 
very  good  agreement  between  the  profiles  has  been 
obtained.  In  the  case  of  highly  strained  flame  the 
temperature  and  major  species  profiles  were 
predicted  well,  while  the  minor  species  were 
overpredicted  in  the  products  region  far  behind  the 
flame.  The  reason  for  these  results  is  the  absence  of 
the  momentum  equation  in  our  model.  Although 
the  momentum  equation  is  redundant  along  the 
stagnation  streamline,  it  provides  means  by  which 
the  boundary  condition  eb  =  £«  V  Pu  Ipb  . 

which  relates  the  strains  on  the  reactants  and 
products  sides,  can  be  imposed.  Our  model  does  not 
satisfy  this  boundary  condition,  accounting  for  the 
density  variation  in  the  reaction  zone  by  the 
combustion  velocity  only,  and  hence  it  is  does  not 
allow  the  strain  to  vary  throughout  the  flame 
structure.  Thus,  effectively,  while  the  flame 
structure  is  resolved  correctly,  the  strain  rate  on  the 
products  side  in  our  model  is  lower  than  its  actual 
value. 

2.  The  transient  dynamics  of  lean  and 
stoichiometric  flames,  in  response  to  a  sudden 
change  in  strain,  have  been  analyzed.  It  has  been 
found  that  the  dynamics  is  determined  by  the  heat 
release  rate  in  the  flame.  Therefore,  the 
stoichiometric  flame  which  has  the  highest  heat 
release  rate  has  the  smallest  response  time.  For 
leaner  mixtures,  flames  respond  more  slowly  and 
their  time  constants  increase. 

3.  For  a  given  equivalence  ratio,  the 
response  time  first  increases  with  strain,  reaching  a 
maximum  at  approximately  700  1/s,  then 
monotonically  decreases.  Two  regimes  of  flame 
response  have  been  identified.  In  the  lower  strain 
regime,  the  flame  is  convected  by  the  flow  while 
its  structure  remains  almost  unchanged  and  the 
flame  dynamics  is  governed  by  the  thermal- 
diffusion  time  scale,  which  increases  with 
equivalence  ratio.  In  the  second  regime  the 
influence  of  the  flow  on  the  flame  structure  is 
dominant,  the  flow  time  scale  determines  the  flame 
dynamics  and  makes  it  independent  from  the 
equivalence  ratio. 

4.  In  turbulent  combustion  simulation, 
our  results  show  that  only  near-stoichiometric 
flames  demonstrates  instantaneous  response  to 


typical  strains  and  can  be  treated  as  steady. 
Otherwise,  the  unsteady  local  flame  response  must 
be  considered  in  detail. 
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Wmix  -  the  molecular  weight  of  mixture, 

Hk  -  the  enthalpy  of  the  k  -th  species,  molar 
units, 

Vkr  the  diffusion  velocity, 

Vcor  -  the  correction  velocity, 

kr,  i  -  the  thermal  diffusion  coefficient  of  the 
k  -th  species, 

(x,  y,  t)  -  physical  coordinates, 

( y.  t)-  coordinates  of  the  first  transformed  domain, 
x)  -  computational  coordinates, 
ef  0  -  the  strain  rate, 

L,,  f  the  Lewis  number  of  the  k  -th  species, 

a-  the  thermal  diffusivity  of  mixture, 

6-  the  non-dimensional  temperature, 
ocki-  the  third-body  coefficient  of  the  k  -th  species 
in  reaction  i, 

Palm-  the  pressure  of  one  attnosphere. 


NOMENCLATURE 
p  -  the  mixture  density, 

V  =  (u,  v)-  the  total  velocity  of  fluid, 

V  -  (  UkV )  -  the  inviscid  irrotational  velocity 

component, 

^c  =  ( Vc the  combustion  generated 
velocity, 

Yk  -  the  mass  fraction  of  the  k  -th  species, 

Xk-  the  mole  fraction  of  the  k  -th  species, 

[Xkj-  the  molar  concentration  of  the  k  -th  species, 
Dk  -the  mixture  averaged  diffusion  coefficient  of 
the  k  -th  species, 

d)k  -  the  production  rate  of  the  k  -th  species, 

Wk  -  the  molecular  weight  of  the  k  -th  species, 

Cp  -  the  constant  pressure  specific  heat  of  mixture, 
mass  units, 

Cp,  k  -  the  constant  pressure  specific  heat  of  the 
k  -th  species,  mass  units, 

X'+i  -the  total  number  of  chemical  species, 

r  -  the  thermodynamic  temperature, 

p  -  the  pressure, 

X  -  the  thermal  conductivity  of  mixture, 

Ru- the  universal  gas  constant. 


SUBSCRIPTS 

u-  unbumt  mixture, 
b-  burnt  mixture, 
scl-  scale  value, 
k-  species  number, 
n-  time  layer  number, 
mix-  mixture. 
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Figure  1.  Flame  propagation  in  a  stagnation  point  flow. 
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Figure  2.  A  comparison  between  the  computed  temperature 
and  mole  fraction  profiles  for  a  stoichiometric,  unstrained 
flame  obtained  by  our  model  (solid  line)  and  by  Smooke's 
model  (dashed  line). 
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Figure  3.  A  comparison  between  the  computed  mole 
fraction  profiles  for  a  stoichiometric,  unstrained  flame 
obtained  by  our  model  (solid  line)  and  by  Smooke's 
model  (dashed  line). 
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Figure  4.  A  comparison  between  the  computed  mole 
fraction  profiles  for  a  stoichiometric,  unstrained  flame 
obtained  by  our  model  (solid  line)  and  by  Smooke's 
model  (dashed  line). 
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Figure  5.  A  comparison  between  the  computed  temperature 
and  mass  flection  profiles  for  a  stoichiometric,  strained 
flame  employing  our  model  (solid  line)  and  the  model  of 
Rogg  (dashed  line);  strain  rate  «  1,000  1/s. 


Figure  6.  A  comparison  between  the  computed  mass 
fiaction  profiles  for  a  stoichiometric,  strained  flame 
employing  our  model  (solid  line)  and  the  model  of 
Rogg  (dashed  line);  strain  rate  »  1,000  1/s. 
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Figure  7,  A  comparison  between  the  computed  mass 
fraction  profiles  for  a  stoichiometric,  strained  flame 
employing  our  model  (solid  line)  and  the  model  of 
Rogg  (dashed  line);  strain  rate  =  1,(XX)  1/s. 


Figure  9.  The  steady-state  heal  release  rate  of  a  premixed, 
mediane-air  flame  as  a  function  of  the  strain  rate  and  the 
equivalence  ratio. 


Figure  11.  The  response  time  of  a  premixed,  methane- 
air  flame  as  a  function  of  the  applied  strain  and  the 
equivalence  ratio. 
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Figure  8.  A  comparison  between  the  computed  mass 
fraction  profiles  for  a  stoichiometric,  strained  flame 
employing  our  model  (solid  line)  and  the  model  of 
Rogg  (dashed  line);  strain  rate  =  \fiO0  1/s. 


Fi^re  10.  A  typical  response  of  a  premixed, 
stoichiometric,  atmospheric  pressure,  methane-air 
flame  to  a  step-wise  change  in  the  strain  rate; 
strain  rate  =  1,000  1/s. 


APPENDIX  m 


A  Uniform  Strain  Model  of  Elemental  Flame  for  Turbulent 


Combustion  Simulation 


Constantin  Petrov  and  Ahmed  Ghoniem, 
Massachusetts  Institute  of  Technology, 
Cambridge  MA  02139,  USA 


I.  Introduction 

The  paper  describee  pan  of  our  effort  to  develop  an  unsteady  strained  laminar  flan*  mode 
with  detailed  chemistry  and  nansport  which  can  be  used  as  a  combustion  submodel  in  a  nnbulent 
combustion  simularion.  The  combustion  submodel  is  based  on  nacking  the  flame  surface  as  tt  is 

expanded,  convoluted  and  strainedby  the  turbulent  structures.  The  flame  surface  is  divided  mm  a 

urge  number  of  elemental  flames,  whose  burning  rate  (diffusion  flame)  or  burning  velocity  ( 
premixed  flame)  is  computed  using  the  elemental  flame  structure  solution.  In  order  m  pmdict 
correcUy  the  burning  tate  and  the  poUutants  fotmation  in  a  turbulent  reacting  flow,  each  flame  and 
iu  interactions  with  other  elemental  flames  must  be  modeled  adequately.  i.e.  all  nnportant 
phenomena  which  affect  instantaneous  flame  structure  must  be  accounted  for  in  numerical  model. 

It  has  been  detemrined  in  alarge  number  of  previous  studies  ( see.  e.g.,  a  review  by  Peters 
(1) )  that  the  strain  rate,  unbalanced  molecular  diffusion  and  curvature  effects  can  significantly  alter 
dm  burning  rate.  In  our  previous  papers  (  Refs.  t2H4)  )  we  established  that  a  flame  does  not 
mspond  instantaneously  m  a  rapid  change  in  dm  snain  rate,  ftequently  making  the  quasr-steady 
assumption  inappropriate.  Similar  results  regarding  the  importance  of  transient  response  were  also 
obtained  in  Reft.  15]  -  [8].  The  requirement  of  predicnng  the  pollutants  formation,  such  as  CX)  and 
NOx.  necessHa^  appBcatkm  of  complex  chemistry  or  reduced  kinetics  mechanisms.  Furthermore, 
at  some  locations  of  a  nonptemixed  combustor  premixed  or  pamally  premixed  combustion  occurs 

and  must  be  simulated  using  the  same  flame  model 

Summarizing,  a  model  of  elemental  flame  must  account  for  the  foUowing  physical 
phoiomena :  (1)  the  strain  ram.  (2)  unsteadiness  of  flame-flow  interaction.  (3)  compl«t  chemistry 
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and  transport  It  must  be  able  to  (4)  predict  diffusion,  premixed  and  partially  premixed  flames  at 
the  same  time.  As  far  as  the  numerical  requirements  are  concerned,  the  model  should  be  (5) 
numerically  effldent  and  (6)  robust  in  order  to  be  able  to  predict  such  extreme  phenomena  as 
ignition  and  extinction.  If  these  requirements  are  to  be  met,  the  influence  of  curvature  on 
combustion  should  be  accounted  on  a  global  scale  by  combining  elemental  flames  and  modeling 
their  interaction.  This  problem  is  beyond  the  scope  of  our  paper. 

The  specifications  of  the  model  are  implemented  as  follows.  The  effect  of  the  strain  rate 
exerted  by  the  flow  on  the  flame  is  accounted  for  by  considering  the  widely  used  configuration  of 
flame  propagating  in  a  stagnation  point  flow.  After  testing  the  shooting  technique  and  finite- 
difference  boundary-value-problem  technique  of  solving  the  equations  governing  flame 
propagation  in  this  type  of  flow,  we  chose  to  solve  the  transient  problem  using  an  implicit  finite 
difference  scheme.  This  approach  satisfies  requirements  (l)-(4)  and  (6).  As  far  as  the  requirement 
(5)  of  numerical  efficiency  is  concerned,  we  found  that  the  rapid  advances  in  the  computer 
technology  and  the  robusmess  of  the  time  integration  procedure  in  the  region  near  flammability 
limits,  ignition  and  in  the  case  of,  for  example,  high  pressure,  outweighs  perceived  high 
computational  cost  as  compared  to  the  finite-difference  boundary-value-problem  techniques  used 
by  Smooke  [10],  Rogg  [1 1],  Kee  et.  al.  [16]. 

In  this  paper  we  propose  two  approaches  to  compute  the  structure  of  an  elemental  flame. 
Both  are  able  to  predict  unsteady  one-dimensional  premixed  as  well  as  diffusion  flames.  In  the 
first,  which  we  call  the  Unsteady  ^.trained  Elemental  flame  (  USEF  ),  we  integrate  the  unsteady 
continuity,  mcxnentum,  energy  and  species  mass  fraction  equations  governing  flame  propagation  in 
a  stagnation  flow.  We  use  this  solution  as  a  benchmark  to  compare  with  the  predictions  of  the 
Uniformly  ^trained  j^emental  flame  (  UniSEF  )  model  in  which  we  assume  a  uniform  effective 
strain  throughout  the  flame  structure.  This  simplification  allows  us  to  decouple  the  energy  and 
mass  fraction  equations  from  the  momentum  equation  and  to  transform  the  former  into  a  set  of 
reaction-  diffusion  equations.  The  advantage  of  solving  the  reaction-  diffusion  equations  governing 
UniSEF  instead  of  the  full  system  governing  USEF  is  the  robustness  of  the  numerical  integration 
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procedure  used  in  the  former  since  the  troublesome  convective  terms  are  eliminated  from  the 
system.  Another  advantage  is  the  reduced  number  of  governing  equations.  While  this  may  have 
little  significance  in  the  case  of  detailed  chemical  kinetics  since  most  of  the  run  time  is  spend 
calculating  the  Jacobian  and  the  chemical  source  terms,  it  could  be  important  if  a  reduced  chemical 
kinetics  mechanism  with  a  smaller  number  of  species  is  adopted. 

The  paper  is  organized  as  follows.  In  Section  II,  we  present  the  formulation  and  the 
reduction  to  the  UniSEF.  In  Section  IE,  we  describe  the  numerical  procedure  used  to  integrate  the 
equations  in  both  cases.  In  Section  IV,  we  compare  the  solution  of  USEF  with  these  of  Ref.  [10], 
and  investigate  the  error  in  UniSEF  by  comparing  its  results  with  those  of  USEF  in  terms  of  the 
steady  state  premixed  and  diffusion  flame  structures  and  their  time  responses. 

n.  Formulation 

In  the  USEF  model,  we  integrate  a  reduced  version  of  the  unsteady  equations  governing 
the  strained  flame  propagation  along  the  stagnation  streamline  of  a  stagnation  point  flow  (12.11- 
14)  of  Ref.  [9].  We  start  with  the  general  unsteady  equations  for  mass,  momentum,  chemical 
species  and  energy : 


dp  ,  d(p  u)  .  d(p\)  _ 
dt  dx  dy  “ 
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where  x  andy  axes  are  directed  parallel  and  popendicular  to  the  flame,  respectively,  r  is  the  time, 
u  and  V  are  the  X  and  y  velocity  components,  p  is  the  density,  p  is  the  pressure,  M  is  die 
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molecular  viscosity,  Yk  is  the  mass  fraction,  Dk  is  the  diffusion  coefficient,  Wk  is  the  molecular 
weight,  Wk  is  the  production  rate,  T  is  the  temperature,  Cp  is  the  specific  heat  of  mixture  in  mass 

units,  A  is  the  thermal  conductivity  of  mixture,  Hk  is  the  total  enthalpy  in  mole  units.  The  system 
of  equations  is  closed  by  the  ideal  gas  law.  We  neglected  the  enthalpy  flux  term  in  the  energy 
equation  and  the  thermal  diffusion  velocity  in  the  species  equations.  It  has  been  shown  by  Smooke 
[10],  and  confirmed  by  us,  that  these  terms  are  relatively  unimportant. 

The  free-stream  velocity  outside  the  boundary  layer  is  given  by 


Uoo-e(t)  X 


(5) 


where  e.  is  the  strain  rate.  In  the  outer  layer  the  y-derivatives  in  the  momentum  equation  arc 
negligible  and  we  can  write 


d  Ugo 


du„ 
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p  Moo 
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pu„e 


which  allows  us  to  express  the  pressure  gradient  as  a  function  of  the  imposed  velocity. 
Considering  the  solution  of  the  above  equations  along  the  stagnation  streamline,  x  =  0,  whCTe 

M  =  0,  —  =  0  and  introducing  the  notations  of  Ref.  [2]  f/  =  «  /  «oo »  V  =  p  v  ,  we  rewrite 
dx 

Eqs.  (1)  -  (4)  as-follows 
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In  the  momentum  equation  there  are  two  terms  containing  the  time  derivative  of  strain.  Since  the 
density  inside  the  flame  is  usually  much  lower  than  the  density  outside  and  U  is  of  order  of  unity, 
the  following  inequality  holds 


p  U  Uoo 


]  de 
e  dt 
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poo  Moo 


2_  de 
^Jt 


In  the  vicinity  of  the  stagnation  streamline  Uoo=  e  x  0  .  It  can  be  demonstrated  using 
dimensional  analysis  that  close  to  the  stagnation  streamline  all  convective  terms  in  the  x-  direction 
are  small  compared  to  the  corresponding  terms  in  the  y-direction.  Dividing  the  momentum  equation 
by  Moo  ,  we  finally  obtain 
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d  U 
dt 
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dt 
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Gose  to  the  stagnation  streamline,  the  remaining  governing  equations  take  the  form 
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In  the  two  jet  configuration,  the  following  boundary  conditions  apply 
at  y=  +  oo  :  Yii=  Yk,+oo,T  =  T+«, ,U  =  1 

at  y  —  -  :  Yi^=  Y^^ .00 ,T  =  T.00 ,  U  =  V p+oJ  P-oo  (14) 

at  y  =  0  :V  =  0 

In  the  Tsuji  configuration,  we  have 

at  y=  +  00  :Yk=  Yk,+c^,T  =  T+00,  U  =  1  ^ 

at  y  —  y^f  .*  Yk  —  Yk^  -  {p  D  dYjjBy  >T  —  ,  t/>v  =  0  ^  (15) 

where  subscript  “w”  indicates  the  values  taken  at  the  burner  surface.  The  mass  flux  at  this  surface, 
is  considered  to  be  a  known  value.  Both  configurations  are  shown  in  Fig.  1.  The  above 
equations  and  the  boundary  conditions  can  be  applied  to  the  strained  diffusion  and  pretnixed 
flames.  Next,  we  formulate  the  alternative  UniSEF  model 

The  UniSEF  model  is  based  on  the  approach  suggested  by  Carrier  et.  al.  [12]  for  diffusion 
flames  with  infinite  rate  chemistry,  extended  by  Rutland  and  Ferziger  [13]  to  the  unity  Lewis 
number  premixed  flame  with  single  step  chemistry,  by  Cetegen  and  Bogue  [7]  to  a  diffusion  flame, 
by  Ghoitiem  etf^  [6]  to  a  diffusion  flame  with  single  step  chemistry  and  by  Petrov  and  Ghoniem 
[2]  to  the  nonuniQf  Lewis  number  premixed  flame  with  single  step  chenustry,  and  to  the  complex 
chemistry  prerruxed  flame  [3].  The  solution  procedure  for  a  generic  strained  (i.e.  premixed  and/  or 
diffusion)  flame  is  identical  to  that  described  in  detail  in  Ref.  [3].  Eqs  (1)*(4)  written  along  the 
stagnation  streamline  take  the  form : 
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where  x-momentum  equation  is  eliminated  because  it  is  trivially  satisfied.  We  decompose  the 
velocity  field  into  two  components : 


V 


=  Vp  +  Vc 


(19) 


where  v  is  the  total  velocity,  Vg  is  the  combustion  generated  velocity,  and  is  the  externally 
imposed  velocity  which  satisfies  the  following  incompressible  continuity  equation: 


(20) 


If  the  externally  imposed  inviscid  iirotational  velocity  field  ,Vp,  is  specified  as 


vp-(e^(t)x.-e^(t)y) 


(21) 


where  ^  is  the  ^ective  strain  rate,  then  the  following  transformations  ( see  Ref.  [3] ),  where 
the  overbar  denotes  the  density  normalized  by  its  value  at  positive  infinity,  p+  oo^ 


{ 


y  =  exp(\  e^(f)  dY)  y  sB(t)y,  t  =  r 


T  =  t 


(22) 
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eliminate  the  Vp  d!  dy  operators  from  the  governing  equations.  The  uniform  strain  model  assumes  a 
constant  value  of  strain  across  the  flame,  ^eff.  However,  in  the  exact  solution,  the  values  of  strain 
change  inside  the  flame  structure  from  emax  and  Emin  which  can  be  significantly  different.  The 
exact  steady  state  profiles  of  temperature,  heat  release  rate  and  normalized  strain  rate  U  in 
premixed  and  diffusion  flame  are  shown  in  Fig.  2.  The  location  of  the  maximum  heat  release  rate  is 
denoted  by  a  vertical  dashed  line.  We  found  that  an  effective  strain  rate  equal  to  the  strain  inside  the 
reaction  zone,  i.e.  the  value  of  strain  at  the  intersection  of  the  dashed  line  and  the  strain  rate  profile 
significandy  improves  the  predictions  of  the  UniSEF  solution.  For  a  premixed  and  diffusion 
flame,  this  value  is  approximately  equal  to  average  of  the  minimum  and  maximum  values  of  strain 
in  the  domam,  Sffidx  and  £fnin*  t.e. 


~  (^max  +  ^min  )! 2  —  £+eo  (  1  (  p+eo  I Pb  )  ^  ^23) 

where  emin-  and  Pb  corresponds  to  the  adiabatic  flame  temperature.  In  the  derivation  of  this 
relationship  the  momentum  equation  we  used  and  the  first  and  second  derivatives  of  U  at  the 

location  of  ^max  assumed  to  be  equal  zero. 

Mapping  the  governing  equations  onto  the  ( t)  domain,  we  get  the  final  form  of  the 
UniSEF  equations; 


^  =:  ^  +  hQ(6kWk  ^ 
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(25) 


These  equations  are  non-dimensionalized  using  a  length  scale  based  on  the  thermal  diffusion 
length,  i.e.  Jscl  -  V  C4  « tjc/.  where  ct).  oo  is  the  thernml  diffusivity  coefficient  of  the  mixture  at 
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infinity,  (,<./  =  10 '^s .  The  non-dimensional  temperature  is  0  =  T  /  7+ 

.  andi  =  r/ tsci,-  The  boundary  conditions  remain  the  same  as  in  Eqs.  (14)  or  (15) 
with  the  exception  of  those  for  U  and  V  which  are  not  required  anymore. 

In  order  to  solve  the  governing  equations  of  USEF  and  UniSEF  a  transport  model  is 
necessary  to  evaluate  die  mixture  conductivity  and  the  diffusion  coefficients  of  chemical  species. 
We  adopted  a  simplified  approach  of  Smooke  [  10] : 

(a)  the  Lewis  numbers  of  different  species  are  constant  in  y  -  direction 


=  const 

f 


(26) 


(b)  the  following  curve  fit  can  be  used  to  calculate  the  ratio  X/  Cp 
XlCp=  2J8  10 -^(71298  )  (CGS  units) 


and 

(c)  the  Prandtl  number  is  constant  and  equal  to  0,75 


By 


(11 


d  U 
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The  validity  of  ftese  assuaq>tions  has  been  extensively  tested  and  confirmed  by  us  using  Chemkin 
Transport  Properties  Package  [14].  However,  Smooke  [10]  demonstrates  that  this  simple  transport 
model  leads  to  a  lower  value  of  the  extinction  strain.  We  exercised  caution  in  applying  this  model 
and  tested  the  results  produced  by  USEF  against  those  obtained  using  the  con^lex  transport 
version  of  the  model  which  is  also  available. 
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The  starting  chemistry  models  of  methane  oxidation  due  to  Smoke  [10]  is  used  in  the 
paper.  It  contains  25  reactions  and  16  species. 

V 

UL.  Numerical  Solution  Method 

The  numerical  solution  of  Eqs  (10)-(13)  and  (24)-(25)  with  boundary  conditions  (14)  or 
(15)  is  obtained  using  an  implicit  Euler  method  and  Newton-like  iterations.  Here  as  an  example  we 
present  the  finite-difference  approximation  of  Eq.  (13).  Other  equations  are  descretized  similarly. 
Assume  that  all  dependent  variables  are  known  at  time  layer  n  .  In  the  implicit  Euler  method,  the 
time  derivatives  are  approximated  as 
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dT 

dt 


'T  n+1  'T  n 

iJ  .  '.// 

At 


(28) 


where  subscripty  indicates  the  grid  point  and  the  superscript  n  denotes  the  n-th  time  layer.  At  is 
the  time  step. 

An  approximation  of  the  convective  term  in  USEF  is  extremely  important  for  the  accuracy 
and  stability  of  the  simulation.  In  our  method,  we  use  a  central  difference  formula 
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where  hj  =  yj  ,  to  calculate  the  first  derivatives  in  the  convective  terms  wherever  it  is 
possible .  In  the  regions  of  potential  numerical  instability  we  use  the  following  upwind  derivatives 
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In  these  regions  the  first  derivative  of  temperature  repeatedly  changes  sign.  The  diffusion  and  the 
source  terms  in  USEF  and  UniSEF  are  approximated  in  the  usual  way 


,  'rn+l  rpn+l 
2  rt+/  I  j  *  ij-l 

hj.i 


(32) 

where  subscript  j+1/2  indicates  the  value  taken  half-way  between  the  j-th  and  the  (j+l)-th  mesh 
points. 

The  resulting  system  of  non-linear  algebraic  equations  is  solved  by  a  damped  Newton 
algorithm.  In  the  USEF  noodel,  we  define  the  i-th  Newton  iteration  of  the  partial  (  without  V  - 

component )  solution  vector  0”"^^  at  the  (n+1 )  time  layer  as  follows 
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In  the  UniSEF  model,  the  solution  vector  has  no  f/  and  V  -  components.  When  an  arbitrary 
solution  vector  is  substituted  into  the  finite-difference  analog  of  Eqs  (10)  -  (13)  and  (24)-(25),  tire 
right-hand-sides  of  these  equations  are  not  equal  to  zero  as  they  would  if  the  true  solution  were 


substituted.  We  denote  this  residual  vector  by  P.  The  purpose  of  the  Newton  iterations  is  to  find  a 

veaor  which  satisfies  the  following  equation 

\ 

F  ( =  0 

with  acceptable  accuracy. 

The  initial  guess  of  the  solution  vector  is  obtained  by  extrapolating  the  converged  solution 
vectors  from  two  previous  time  steps.  Once  the  initial  guess  is  obtained,  the  next  iteration  for  the 
solution  vector  at  the  (n+1)  time  layer  is  calculated  using  a  damped  Newton-  like  algorithm 

= ■x<‘> 

W*/  ,  (33) 

where  BF  /  is  the  numerically  evaluated  Jacobian  matrix , 

(BF  (  F(  +  S)  -  5  ,  S  =  0.0001 

and  X  is  a  damping  coefficient  of  order  of  unity.  For  the  given  system  of  equations  the 
Jacobian  matrix  has  a  block  -  tiidiagonal  smicture.  The  system  of  linear  algebraic  Eqs  (33)  is 
solved  numerieii^  using  a  standard  numerical  procedure  (  see  Ref.  [15]  ).  The  Newton  -  like 

Ji)\  ^n-¥l,  (i+1)  .n+1.  (i)  I 

iterations  ctmtinue  until  a  norm  of  the  conection  vector ,  =  y  ~  V  I  is  smallCT 

than  a  user  specified  tolerance  parameter.  The  norm  is  calculated  as  follows 
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For  every  time  step,  a  single  iteration  was  often  sufficient  for  convergence.  In  USEF  ,  in 
the  case  of  two  jet  configuration,  once  the  profiles  of  chemical  species,  temperature  and  U-  velocity 
are  obtained  at  the  ('«+/)  time  layer,  the  mass  flux  is  calculated  by  integrating  the  continuity 
equation 


dt 
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pUe  + 


dV 

dy 


from  the  stagnation  point  y  -  0  ,  where  V  =  0,  to  the  left-  and  right-  hand  boundaries  of  the 
computational  domain.  The  dV !  dy  derivative  is  approximated  by  an  upwind  scheme 


V///  -  Vf 


pUe)f:/a  =0,y>0 


pUe)]^l}2  =0,  y  ^0 


13 


In  the  Tsuji  flame  configuration  the  continuity  equation  is  integrated  from  the  burner  surface,  y^, 
where  the  value  of  is  known,  to  infinity.  In  UniSEF  ,  only  the  temperature  and  species 
equations  are  solved  and  no  integration  of  the  continuity  equation  is  required. 

After  the  complete  solution  vector  is  obtained  in  both  models,  the  computational  mesh  is 
adapted  following  a  procedure  described  in  detail  in  Ref.  [16].  At  every  time  step  the  thermal 
thickness  of  the  flame  is  compared  with  the  distance  from  the  flame  to  the  upstream  and 
downstream  boundaries  of  the  computational  domain.  If  either  of  these  distances  is  smaller  than 
two  flame  thickness,  the  computational  domain  is  regridded  by  moving  the  boundaries  away  from 
the  flame.  The  time  step  of  the  calculations  is  taken  to  be  the  minimum  of  the  characteristic  time 
scales  of  species  production  and  is  bounded  to  be  lower  than  l.Oe-06  s. 


rv  Results  and  Discussion 

In  Fig.  3  the  steady  state  structures  of  an  atmospheric  pressure,  stoichiometric,  diffusion, 
methane-air  flame  with  a  strain  rate  e+oc  =  100  1/s  in  the  Tsuji  configuration  obtained  using  USEF 
(solid  lines  )  and  a  phase-space,  pseudo-arc-length  continuation  method  [10]  (dashed  lines)  are 
compared.  The  same  chemical  kinetics  and  transport  models  are  used  in  both  simulations.  The 
profiles  of  species  and  temperature  are  shown  as  a  function  of  the  mixture  fraction,  Z,  defined  as  a 
solution  of  the  mixture  fraction  transport  equation 

"  dy  Sy\cp3y]  ,  (34) 

where  the  mixture  fraction  is  equal  to  unity  on  the  oxidizer  side  and  to  zero  on  the  fuel  side.  We 
observe  an  excellent  agreement  in  terms  of  the  shapes,  locations  and  peak  values  of  the  temperature 
and  chemical  species.  Some  discrepancies  in  the  CH4  profile  around  Z  =  0.07  are  enhanced  by 
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the  choice  of  the  mixture  fraction  as  an  independent  variable  and,  probably,  due  to  a  high  numerical 
diffusion  of  the  upwind  discretization  scheme  applied  in  some  grid  points  of  our  computational 
domain.  ! 

The  above  test  run  serve  to  demonstrate  that  the  USEF  model  is  capable  of  predicting 
accurately  the  steady  state  structure  of  strained  flame  and  that  it  can  be  used  as  a  benchmark 
solution  for  a  comparison  with  the  simplified,  uniform  strain  UniSEF  model. 

Comparison  of  USEF  and  UniSEF 

The  USEF  and  UniSEF  models  are  compared  in  terms  of  the  transient  dynamics  and 
Steady  state  premixed  and  diffusion  flame  structures.  In  this  paper  an  atmospheric  pressure, 
premixed  flame  with  a  strain  rate  £+  oe  =  1000  1/s  and  an  atmospheric  pressure,  diffusion  flame 
with  a  strain  rate  =  200  1/s  are  used  as  examples.  In  both  test  cases,  the  flames  are  simnlateH 
in  the  two  jet  configuration  and  their  steady  state  structures  are  shown  in  Figs.  7  and  8, 
respectively. 

Before  we  compare  the  transient  dynamics  and  steady  state  structures  of  USEF  and 
UniSEF,  it  is  worthwhile  to  understand  how  the  assumption  of  uniform  strain  affects  the  velocity 
and  mass  flux  fields  in  both  flames.  Fig.  4  shows  steady  state  heat  release,  mass  flux  and  y- 
velocity  profiles  obtained  using  both  solutions  in  the  premixed  flame  case.  Fig.  5  shows  the  same 
profiles  in  the  diffusion  flame.  In  both  figures  all  UniSEF  profiles  are  drawn  by  dashed  lines, 
USEF  profiles  -  by  solid  lines.  We  note  that  near  the  peak  values  of  the  heat  release  rate  profiles 
the  mass  flux  and  y-vclocity  profiles  of  the  exact  and  simplified  solution  are  very  close.  This  is  due 
to  the  choice  of  Ae  unifonn  strain  (  see  Eq.  (23)  )  as  the  exact  strain  inside  the  reaction  zone. 
Q>mbustion  inside  the  reaction  zone  in  UniSEF  occurs  at  almost  the  same  conditions  as  in  USEF. 
In  the  next  section  we  will  see  that,  althou^  outside  the  reaction  zone  the  mass  flux  and  velocity 
profiles  of  the  two  solutions  differ  significantly,  their  closeness  inside  the  zone  leads  to  a  very 
good  agreement  in  terms  of  the  transient  dynamics  and  steady  state  structures. 
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Premixed  Flame 


We  start  (fee  comparison  of  the  models  with  the  premixed  flame.  Both  solutions  are  initiated 
with  identical  arbitrary  profiles  of  the  dependent  variables  constructed  from  the  hyperbolic  tangent 
and  bell-shaped  exponential  functions  or  a  combination  of  both.  The  values  of  the  variables  far 
behind  the  flame  on  the  products  are  set  equal  to  the  equilibrium  values  for  a  given  pressure  and 
chemical  kinetics  mechanism.  The  peak  values  of  radical  profiles  are  of  the  order  of  typical  values 
encountered  in  premixed  methane-air  flames.  The  initial  thickness  of  all  profiles  is  6  mm,  the  initial 
location  of  flame,  defined  as  the  location  of  the  maximum  fuel  consumption,  is  1  mm  from  the 
stagnation  point.  The  chemical  kinetics  and  transport  models  of  Ref.  [10]  are  used  in  all 
simulations.  At  time  t  =  0  s.  the  initial  profiles  are  introduced  into  the  models  and  the  strain  rate 
£4.  ee  is  instantaneously  put  equal  to  1000  1/s. 

In  Fig.  6a  the  USEF  time  histories  of  the  burning  velocity  and  flame  location  are  shown  by 
lines  with  empty  symbols,  the  UniSEF  time  histories  -  by  lines  with  solid  symbols.  Fig.  6a 
demonstrates  virtually  identical  time  response  of  USEF  and  UniSEF  to  a  step-wise  change  in  the 
strain  rate.  This  figure  clearly  show  that  the  transient  response  of  a  premixed  flame  can  be  modeled 
using  UniSEF  .  In  the  case  considered,  the  methane-air  premixed  flame  response  time  of 

0  ( 10'^ )  s.  This  could  be  several  times  higher  than  a  typical  time  step  in  a  turbulent  flow  making 
the  quasi  -  steady  assumption  invalid.  A  comparison  of  the  quasi-steady  approximation  and  die 
transient  approach  is  made  later  in  the  paper  using  the  diffusion  flame  case  as  an  example. 

The  steady  state  profiles  of  temperature  and  major  species  produced  by  USEF  (solid  line) 
and  UniSEF  Ctlashcd  Une)  for  the  same  flame  are  compared  in  Fig  7.  The  shapes  of  the 
temperature,  major  and  minor  species  profiles  are  virtually  identical,  although  the  location  of  the 
UniSEF  flame  is  further  from  the  stagnation  point  than  that  of  the  USEF  solution  due  to  a  much 
lower  absolute  value  of  the  uniform  strain  on  the  products  side  which  dominates  over  a  somewhat 
higher  value  of  the  strain  on  the  reactants  side  ( see  Fig.  2a) .  Peak  values  of  radicals  also  agree 
well. 
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We  conclude  that  in  the  case  of  the  premixed  flame  the  uniform  strain  UniSEF  model 
produces  the  steady  state  flame  structure  and  time  response  which  agree  well  with  the  results  of  the 
exact  USEF  model. 

Dijfusion  Flame 

The  same  two  approaches  are  used  to  smdy  the  structure  and  the  time  response  of  a 
diffusion  flame  subjected  to  a  strain  e+  «.  of  200  1/s.  Again,  both  simulations  start  with  the  same 
set  of  arbitrary  bell-shaped  profiles  of  the  thickness  of  6  mm  and  the  peak  locations  of  1  mm.  The 
ignition  of  the  diffusion  flame  is  achieved  by  introducing  peak  H,  OH  and  O  concentrations 
which  are  an  order  of  magnitude  higher  than  the  typical  values  encountered  in  methane-air 
diffusion  flames.  The  time  responses  are  compared  in  Fig.  6b.  The  burning  rate  in  the  case  of  a 
diffusion  flame  is  normalized  in  such  a  way  that  it  has  the  dimensions  of  velociQr  in  order  to  use 
the  same  characteristic  of  combustion  rate  in  both  premixed  and  diffusion  flame  cases: 


^uel  ^fuel  (fy 


|p+  oo  (  ^fuel,  +  00  -  ^fml,  -oo 


) 


The  response  times,  steady  state  flame  locations  and  burning  velocities  of  the  solutions  agree  quite 

well.  The  response  time  of  the  diffusion  flame  is  again  0(10'^)  s. 

In  Fig.  9  o(»ztparison  of  the  quasi-steady  and  the  transient  approaches  is  shown  in  terms  of 
the  burning  velocity.  The  applied  strain  as  a  function  of  time  is  also  shown,  representing  four 

global  flow  simulation  time  steps  of  10'^  seconds  each.  The  initial  strain  of  100  1/s  is  applied  at 

time  0  which,  after  a  time  step  of  10*^  second,  is  changed  to  200,  400  and,  finally,  to  50  1/s.  At 
time  t  =  0  the  steady  state  profiles  corresponding  to  a  strain  of  100  1/s  are  introduced  into  USEF  . 
During  the  initial  flow  time  step,  when  the  strain  is  100  1/s,  the  burning  velocity  remains 
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unchanged  At  time  2  10'^  s.,  the  external  strain  is  changed  and  the  flame  structure  starts  to  evolve. 
Since  the  flame  response  time  is  much  larger  than  the  flow  time  scale,  the  burning  velocity  does  not 
have  enough  time  to  reach  a  steady  state  value  of  1.02  cm/  s  corresponding  to  the  strmn  of  200  1/s. 
At  the  time  for  the  next  strain  change,  flame  is  still  unsteady,  lagging  behind  the  quasi-steady 
value.  In  Fig.  8  we  demonstrate  that  the  quasi-steady  assumption  can  produce  significant  errors  in 
the  burning  rate  predictions. 

A  comparison  of  the  steady-state  profiles  of  the  temperature  and  species  mass  fractions  is 
shown  in  Fig.  8.  In  the  case  of  a  diffusion  flame  the  overall  flame  structure  produced  by  UniSEF 
agrees  well  with  that  of  USEF. 

In  a  numerical  simulation  of  complex  chemistry  most  of  the  computation  time  is  spend  on 
the  evaluation  of  the  chemical  source  terms.  Elimination  of  the  momentum  and  continuity  equations 
from  UniSEF  is  of  little  value  until  the  total  number  of  chemical  species  is  substantially  reduced. 
Peters  and  Kee  [17]  proposed  a  systematically  reduced  four-step  mechanism  which  is  able  to 
provide  the  essential  properties  of  methane-air  diffusion  and  premixed  flames.  This  mechanism 
makes  use  of  7  so  called  “non-steady”  species,  namely  CH4,  02,  C02,  CO,  H2,  H20,  and  H 
instead  of  the  total  of  16  species  of  the  “starting”  mechanism.  The  other  species  are  assumed  to  be 
in  equilibrium  and  various  algebraic  equilibrium  equations  are  used  to  evaluate  their  concentrations 
in  terms  of  the  concentrations  of  “non-steady”  species. 

Incorporation  of  the  reduced  chemical  kinetics  into  UniSEF  can  potentially  provide  a 
numerically  efficient  and  robust  model  of  elemental  flame  for  turbulent  combustion  simulation. 
However,  several  questions  remain  to  be  answered.  First,  how  the  rather  crude  UniSEF 
assumption  of  uniform  strain  across  the  flame  affects  the  S-curve  and,  in  particular,  the  value  of 
the  extincticm  strain.  Second,  how  close  are  the  transient  responses  of  the  exact  solution  and 
UniSEF  with  reduced  kinetics  to  a  rapid  change  in  strain.  In  the  following  we  will  answer  the  first 
question  and  make  some  comments  on  the  second. 

In  Fig.  10  the  traditional  S-curve  for  an  atmospheric  pressure  methane-air  diffusion  flame 
obtained  using  UniSEF  with  reduced  chemistry  (solid  line  ),  exact  solution  with  full  chemistry 
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[17]  ( dashed  line  with  squares ),  and  exact  solution  with  reduced  chemistry  [17]  (dashed  line  with 
circles)  is  shown.  The  overall  agreement  between  the  exact  solutions  and  UniSEF  is  satisfactory  in 
the  whole  range  strains.  The  extinction  strains  for  UniSEF  is  around  3(X)  1/s,  for  exact  solution  - 
352  1/s,  and  for  exact  solution  with  reduced  chemistry  -  400  1/s.  Thus,  UniSEF  with  reduced 
chemistry  provides  extinction  strain  values  which  are  as  different  from  the  exact  values  as  the 
values  of  any  solution  employing  reduced  chemistry. 

As  far  as  the  unsteady  responses  of  USEF  and  UniSEF  with  reduced  chemistry  is 
concerned,  our  preliminary  results  indicate  that  they  are  very  close  to  the  corresponding  responses 
of  the  models  with  the  “starting”  mechanism.  This  is  rather  surprising  since,  as  it  is  mentioned  by 
Peters  in  Ref.  [18],  in  the  highly  unsteady  initial  stages  of  premixed  and  diffusion  combustion  the 
steady  state  assumptions  of  the  reduction  strategy  are  grossly  violated.  We  plan  to  address,  this 
question  in  the  nearest  future. 

Conclusions 

Two  approaches  were  developed  to  study  the  structure  and  dynamic  characteristics  of 
complex  chemistry,  strained,  premixed  and  diffusion  flames  and  for  the  purpose  of  using  them  as 
elemental  flame  submodels  in  turbulent  combustion  simulation.  In  the  exact  flame  solution  the 
unsteady  temperature,  mass  fractions,  continuity  and  momentum  equations  are  solved  using  an 
implicit  time  -  integration  scheme.  In  the  simplified,  uniform  strain  model  the  governing  equations 
are  reduced  to  the  reaction-diffusion  equations  potentially  providing  a  more  robust  solution 
procedure  and  reducing  the  number  of  governing  equations  by  two.  Since  the  actual  strain  varies 
inside  the  flame  from  Bmax  to  the  average  of  these  values  is  proposed  to  serve  as  the  uniform 
strain  m  the  modeL  For  a  premixed  and  diffusion  flame,  this  assumption  provides  a  remarkable 
agreement  in  the  time  responses  and  the  steady  state  flame  structures  of  the  uniform  strain  and 
exact  solutions.  Incorporation  of  the  Peters  reduced  chemistry  mechaiusm  [17]  into  the  model 
allows  to  take  the  full  advantage  of  the  reduced  number  of  UniSEF  governing  equations  and  to 
obtain  the  extinction  strain  and  S-curve  which  are  satisfactorily  close  to  the  exact  values. 
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Using  the  response  of  the  diffusion  flame  to  strain  as  an  example,  we  demonstrate  that  the 
application  of  the  quasi-steady  assumption,  when  the  instantaneous  value  of  the  burning  velocity 
correspond  to  the^instantaneous  value  of  local  strain,  can  potentially  lead  to  significant  errors  since 
the  flame  response  time  might  be  larger  than  the  flow  calculation  time  step  which  determines  the 
time  interval  between  the  strain  changes. 
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Figure  1.  Two  jet  premixed/  diffusion  and  Tsuji  diffusion  flame  configurations  in  a  stagnation 
point  flow. 

Figure  2.  The  choice  of  effective  uniform  strain  as  the  exact  strain  inside  the  reaction  zone.  USEF 
steady  state  temperature,  strain  and  heat  release  rate  profiles  in  the  premixed  (  strain  =  1000  1/s) 
and  diffusion  (  strain  =  100  1/s)  stoichiometric,  atmospheric  pressure  flames. 

Figure  3.  A  comparison  of  steady  state  (a)  temperature,  major  and  (b)  minor  species  profiles  for  an 
atmospheric  pressure,  diffusion,  stoichiometric,  methane-air  flame  using  the  USEF  model  (solid 
lines )  and  aiclength  continuation  method  (Smooke,  Ref.  [10],  dashed  line);  strain  rate  =  100  1/s, 
chemical  kinetics  mechanism  of  Ref.  [10],  Tsuji  flame  configuration. 

Figure  4.  A  comparison  of  (a)  steady  state  heat  release  and  mass  flux  profiles  and  (b)  steady  state 
heat  release  and  y-velocity  profiles  for  an  atmospheric  pressure,  premixed,  stoichiometric, 
methane-air  flame  using  the  USEF  model  (solid  lines  )  and  the  UniSEF  method  (dashed  line); 
strain  rate  =  1000  1/s,  chemical  kinetics  mechanism  of  Ref.  [10],  two  jet  flame  configuration. 

Figure  5.  A  comparison  of  (a)  steady  state  heat  release  and  mass  flux  profiles  and  (b)  steady  state 
heat  release  and  y-velocity  profiles  for  an  atmospheric  pressure,  diffusion,  stoichiometric, 
methane-air  flame  using  the  USEF  model  (solid  lines  )  and  the  UniSEF  method  (dashed  line); 
strain  rate  =  200 1/s,  chemical  kinetics  mechanism  of  Ref  [10],  two  jet  flame  configuration. 

Figure  6.  A  cco^arison  of  time  histories  of  (a)  the  burning  velocity  employing  USEF  (empty 
symbols)  and  UniSEF  (solid  symbols)  models  for  an  atmospheric  pressure,  premixed, 
stoichiometric,  methane-air  flame;  strain  rate  =  1000  1/s,  and  (b)  the  burning  velocity  employing 
USEF  (empty  symbols)  and  UniSEF  (solid  symbols)  models  for  an  atmospheric  pressure, 
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diffusion,  stoichiometric,  methane-air  flame;  strain  rate  =  200  1/s;  chemical  kinetics  mechanism  of 
Ref.  [10],  two  jet  flame  configuration. 

Figure  7.  Calculated  (a)  temperature,  major  and  (b)  minor  species  profiles  for  an  atmospheric 
pressure,  premixed,  stoichiometric,  methane-air  flame  using  the  uniform  strain  UniSEF  model 
(dashed)  and  the  exact  USEF  model  (solid);  strain  rate=  1000  1/s,  chemical  kinetics  mechanism  of 
Ref.  [10],  two  jet  flame  configuration. 

Figure  8.  Calculated  steady  state  (a)  temperature,  major  and  (b)  minor  species  profiles  for  an 
atmospheric  pressure,  diffusion,  stoichiometric,  methane-air  flame  using  the  uniform  strain 
UniSEF  model  (  dashed  lines  )  and  the  exact  USEF  solution  ( solid  line  );  strain  rate=  200  1/s, 
chemical  kinetics  mechanism  of  Ref.  [10],  two  jet  flame  configuration. 

Figure  9.  A  comparison  of  the  transient  solution  obtained  using  USEF  and  the  quasi-steady 
solution  for  an  atmospheric  pressure,  diffusion,  stoichiometric,  methane-air  flame;  chemical 
kinetics  mechanism  of  Ref.  [10],  two  jet  flame  configuration;  the  strain  history  given  at  the  bottom 
of  the  figure. 

Figure  10.  Calculated  maximum  steady  state  flame  temperature  as  a  function  of  the  inverse  of  strain 
for  an  atmospheric  pressure,  diffusion,  stoichiometric,  methane-air  flame  using  the  uniform  strain 
UniSEF  model  with  reduced  chemistry  (  solid  line  ),  the  exact  solution  with  starting  chemistry 
[17]  (dashed  line,  squares  ),  and  the  exact  solution  with  reduced  chemistry  [17]  (  dashed  line, 
circles );  Tsuji  flame  configuration. 
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While  in  a  number  of  papers  the  reduced  chemistry  approach  is  proven  to  produce  good 
results  for  the  cases  of  steady  state  diffusion  and  unstrained  premixed  flames,  we  are  unaware  of 
any  piblicadon  where  the  method  is  applied  to  study  the  unsteady  strained  premixed  and  diffusion 
flames.  This  is  very  unfortunate  since  the  application  of  it  can  be  potentially  very  helpful  in  the 
simualtion  of  elementary  flame  in  turbulent  combustion.  The  reason  for  the  absence  of  the  unsteady 
applications  of  the  method  is  that,  strictly  speaking,  in  the  highly  unsteady  initial  stages  of 
prentixed  and  diffusion  combustion  its  steady  state  assumptions  are  violated.  Different  chemical 
reactions  are  important  at  that  time.  This  is  mentioned  by  Peters  in  Ref.  [18]. 
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ABSTRACT 

The  transport  element  method  is  described  and  implemented  in  the  simulation  of  the  non- 
premixed  reacting  shear  layer.  The  method,  a  natural  extention  of  the  vortex  element 
method,  resolves  the  variable  density  flow  at  low  Mach  number  and  the  exothermic  reacting 
field.  ITie  effect  of  combustion  on  the  flow  is  accommodated  by  incorporating  a  volumetric 
expansion  velocity  component  and  by  modif5ting  the  integration  of  the  vorticity  equation  to 
include  expansion-related  and  baroclinic  terms.  The  reacting  field  equations  describing  a 
single  step,  irreversible,  chemical  reaction,  are  simplified  by  the  introduction  of  Shvab- 
Zeldovich  (SZ)  conserved  scalars  whose  transport  is  sufficient  to  compute  the  evolution  of 
combustion  in  the  case  of  infinite  reaction  rate.  In  the  case  of  finite  rate  chemistry  the 
evolution  of  one  primitive  scalar,  the  product  mass-fraction,  is  also  computed.  The 
vorticity,  conserved  scalar  gradient  and  product  mass  fraction  are  discretized  amongst  fields 
of  transport  elements.  Their  time  evolution  is  implemented  by  advecting  the  elements  at  the 
local  velocity  while  simultaneously  integrating  their  transport  equations  along  particle 
trajectories.  The  integration  of  the  vorticity  and  the  conserved  scalar  gradient  equations  is 
simplified  using  ideas  fi’om  kinematics.  A  novel  core  expansion  scheme  that  avoids  the 
problems  associated  with  the  conventional  implementation  is  used  to  simulate  diffusion. 
Field  quantities  are  obtained  using  convolutions  over  the  elements.  Results  indicate  that  the 
method  is  able  to  accurately  reproduce  the  essential  features  of  the  flow.  Convergence  of 
the  solution  in  time  is  approximately  linear.  Moreover,  the  finite  reaction  rate  solution  at 
low  Karlovitz  number  bear  strong  similarities  to  that  of  the  infinite  reaction  rate  model. 
This  similarity  is  exploited  in  validating  the  part  of  the  numerical  methodology  related  to  the 
integration  of  the  product  mass-fraction  equation. 


*  Current  address;  Department  of  Mechanical  Engineering,  University  of  Connecticut,  Storrs,  CT06269- 
3139 
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I.  INTRODUCTION 

Generalization  of  the  vortex  element  method  (VEM)  [1]  by  incorporating  the 
physics  of  scalar  transport,  mixing  and  reaction;  what  we  have  called  the  transport  element 
method  (TEM),  has  been  an  ongoing  effort  since  its  inception  [2-9].  The  essential  features 
of  this  process;  the  discretization  of  primitive  variable  gradients,  their  Lagrangian  transport 
and  use  of  flow  kinematics  to  monitor  their  evolution,  offer  significant  advantages  in  the 
simulation  of  reacting  flows.  The  simulation  of  reacting  flows  involve  substantial 
complexities.  These  include  the  presence  and  continuous  mixing  of  many  species,  the 
evolution  of  temporal  and  spatial  scales  which  cover  several  decades,  the  presence  of 
strongly  non-linear  and  coupled  phenomena  and  the  coupling  of  the  combustion  energy 
release  and  the  associated  changes  in  the  density  and  transport  properties  and  the  flow.  The 
latter  requires  the  modification  of  VEM  which  in  its  basic  form  deals  with  incompressible 
uniform  density  flow.  Moreover  reacting  flows  are  characterized  by  thin  reaction  fronts, 
i.e.  regions  of  sharp  scalar  gradients,  which  evolve  in  a  highly  convoluted  manner. 

Early  attempts  to  simulate  reacting  flows  using  a  VEM-based  approach  simplified 
the  reacting  field  through  the  use  of  flame  sheet  models.  In  these  [2]  and  [3],  the  sharpness 
of  the  flame  front  was  exploited  to  represent  it  as  a  sheet  which  encompasses  the  reaction 
and  diffusion  regions  of  the  flame  and  separates  unmixed  scalars.  The  effects  of 
combustion  on  the  flow  were  restricted  to  those  associated  with  the  volumetric  expansion 
by  distributing  volumetric  sources  along  the  flame  sheet. 

To  move  beyond  these  simplified  models,  one  needs  a  detailed  description  of  the 
mixing  field  as  a  first  step  in  the  accurate  representation  of  the  reacting  field.  This  was 
accomplished  by  the  development  of  the  Transport  Element  Method  (TEM)  in  two  [4]  and 
three  dimensions  [5-7].  The  TEM  bears  substantial  similarities  to  the  VEM  to  which  it  is 
fully  compatible.  Essentially,  it  discretizes  the  scalar-gradients  and  evolves  them  by  using 
information  from  the  kinematics  of  material  lines.  The  scalars  are  obtained  from  the 
gradients  via  convolutions  over  their  fields.  A  version  of  the  TEM  incorporating  reaction  in 
a  two-dimensional  temporally  evolving  flow  was  also  proposed  [7]^  Its  implementation, 
however,  in  the  spatially  developing  model  of  the  flow  is  not,  in  general,  possible.  The 
reason  for  this  lies  in  the  temporal  model,  which  due  to  its  utilization  of  periodic  boundary 
conditions  represents  a  mathematical  idealization  of  the  flow  it  attempts  to  represent^  .  The 

1  A  three  dimensional  scheme  was  also  proposed  [8,9]  but  one  in  which  the  scalar-gradient  approach  was 
abandoned  in  favor  of  a  primitive  scalar  one  and  temperature  independent  kinetics  was  assumed. 

2  The  limitations  of  the  temporal  model  in  the  simulation  of  reacting  flows  exceed  those  for  the  non¬ 
reacting  flow  which  have  been  pointed  out  previously  [10].  By  constuction,  the  temporal  model  is  valid 
only  for  very  small  and  very  large  values  of  the  Karlovitz  number.  The  first,  represents  very  fast 
combustion  characterized  by  thin  continuous  flames,  while  for  the  second,  gaseous  fuel  combustion  is 
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TEM  as  presented  in  Ref.  [7]  takes  full  advantage  of  the  simplifications  offered  by  the 
temporal  model  and,  as  a  result,  cannot  be  applied  in  spatially  developing  flows. 

In  this  work,  we  present  a  new  version  of  TEM  which  is  applicable  to  spatially 
developing  flows.  Through  the  introduction  of  SZ  variables,  the  method  is  able  to  use  a 
combination  of  scalar-gradients  and  a  primitive  scalar  to  completely  describe  the  reacting 
field.  Combustion  energy  release  effects  on  the  flow,  at  low  Mach  number,  are  captured 
using  the  variable  density  version  of  the  VEM. 


II.  FORMULATION 

A  two-dimensional  flow  is  considered.  Compressibility  effects  are  restricted  to  the 
low  Mach-number  limit.  Combustion  follows  a  single-step  reaction  which  consumes  two 
reactants,  one  from  each  stream,  to  form  a  single  product  according  to  a  singe  step 
Arrhenius  mechanism.  All  species  behave  as  perfect  gases  with  equal  molecular  weights 
and  equal  and  constant  specific  heats  and  mass  diffusivities.  The  Lewis  number  is  equal  to 
unity  and  molecular  diffusion  is  relatively  small,  i.e.,  the  flow  is  at  high  Reynolds  and 
Peclet  numbers^  .  Accordingly,  the  non-dimensionalized  governing  equations  are: 


dt 
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unsustainable.  Irrespective  of  the  numerical  method  used,  the  model's  implementation  to  flows  of  moderate 
Karlovitz  number  in  which  interaction  of  the  flow  and  reaction  time  scales  is  experienced  and  sharply 
transient  phenomena  such  as  quenhing,  and  re-ignition  take  place  is  problematic. 

^  Most  of  the  assumptions  related  to  the  reaction  model  are  employed  to  reduce  the  computational  effort  and 
are  not  imposed  by  a  fundamental  limitation  in  the  numerical  scheme  (Section  III).  An  exception  to  this  is 
the  equal  diffussivities  assumption  which  enables  the  use  of  the  Shvab-Zeldovich  formulation.  The 
diffussivities,  however,  do  not  necessarily  have  to  be  constant  Moreover,  a  larger  number  of  species  may 
be  considoed  and  the  molecular  masses  may  be  unequal  as  long  as  it  is  assumed  that  the  diffusivides  are  not 
similarly  affected.  This  approach  was  followed  in  the  implementation  of  the  methodology  presented  herein 
for  the  simulation  of  the  axisymmetric  buoyant  plume  [1 1]. 
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♦  ♦  T 

while  p  T  =  1;  Tif  +  ->  (l+(l)*)Tip^  +  Y„  +  Yp  =  1  ^nd  w  =  Yj  Y^  exp(-  . 

In  the  above:  u  =  (u,v)  is  the  velocity  vector  in  a  right-handed  Cartesian  coordinate  system 

+  u  •  V  the  material 

derivative,  p  is  the  pressure,  p  and  T  are  the  fluid  mixture  density  and  temperature 

respectively,  subscripts  f,o  and  p  indicate  the  fuel,  oxidizer  and  product,  respectively,  Tii 

denotes  species  i,  i  =  f,o,p,  Yi  is  the  mass-fraction  for  species  i,  w  is  the  reaction  rate.  Re 

and  Pe  are  the  Reynolds  and  Peclet  numbers,  respectively,  Af,  Qo  and  Ta  are  the  non- 

dimensionalized  frequency  factor,  enthalpy  of  reaction,  and  activation  temperature, 
* 

respectively,  and  (|)  and  <|)  are  the  molar  and  mass  stoichiometry  ratios,  respectively. 

The  equations  are  transformed  into  vorticity-gradient  form,  and  the  velocity  is 
decomposed  into  a  vorticity-induced  solenoidal  component,  u©,  and  two  iirotational 
components;  Ue  due  to  volumetric  expansion, ,  and  Ub  due  to  the  boundary  conditions,: 

u  =  U()i)+  Ub  +  Ue  ,  (4) 

U(o  is  obtained  from  the  definition  of  the  vorticity  and  by  using  a  streamfunction,  \j/,  i.e. 

vV  =  -  £0 ,  Uco  =  Vx(\|/k) ,  (5) 

where  k  is  the  unit  vector  normal  to  the  plain  of  motion.  The  irrotational  components  are 
obtained  from  the  continuity  equation  by  introducing  two  velocity  potentials,  O,  such  that 

V^<De  =  -l^,  Ue  =  V(De  (6) 

p  dt 

and  V^Ob  =  0,  Ub  =  V(Db.  (7) 

The  evolution  of  the  vorticity  is  described  by  taking  the  curl  of  Eq.(2): 

+  (V  u)  cok  =  V^cok  .  (8) 

dt  p2  Re 

For  the  scalar  transport  equations,  Shvab-Zeldovich  (SZ)  conserved  variables,  X 
and  y,  are  introduced  such  that: 

>.  =  Yo-<l)Yf,  Y  =  T--^Yp.  (9) 

1+<|) 


(  d  d  \  H 

X  =  (x,y),  t  is  time,  V  =  gradient  operator,  ^  = 


at 
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These  variables  are  constructed  such  that  their  governing  equations  are  devoid  of  reaction 
source  terms.  If  convection-diffusion  equations  have  the  same  diffusivities  and  boundary 
conditions,  a  reduction  of  SZ  equation  set  to  a  single  equation  can  be  achieved  via 
appropriate  scaling: 

.where =  (10) 

and  X,  1 ,  X.2,  Yi ,  Y2  are  the  scaling  values'* 

In  analogy  to  the  treatment  of  the  flow  equations,  the  S-Z  scalar  transport  equation 
is  recast  into  gradient  form:  g  =  V  P 

^+g  Vu  +  gx(o)k)=^V^g  ,  (11) 

When  the  reaction  rate  is  finite  reactants  do  not  coexist  and  the  reaction  zone 
collapses  onto  a  line.  Under  these  conditions  the  distributions  of  the  S-Z  variables  provide 
a  complete  description  of  the  reacting  field  without  the  need  for  the  solution  of  a  reacting 
scalar.  The  reaction  interface  is  trivially  established  as  the  locus  of  X,  =  0,  while 

\>0  :  Yf  =  0,  Y^  =  X,  Y  =1-X,and  X<0  :  Y,=  -i  Y  =0.  Y  =  1 +i- 

•f  <l> 

When  the  reaction  rate  is  infinite,  the  complete  description  of  the  scalar  field  requires  the 
solution  of  one  reacting  scalar,  e.g.,  the  product  mass-fraction  and 
y  -X.  +  (1-Yp)  X.-Kl)(l-Yp)  Qo  Y 

1+(|)  1-H|)  1  +  (|) 


III.  NUMERICAL  SCHEME 

Discretization:  The  vorticity,  expansion  source,  -  ^  ^ 


P -gradient  and  product-mass- 


fraction  are  represented  by  a  generic  discretization  function.  Thus,  a  property  ^  is 
discretized  among  a  number  of  elements  characterized  by  a  finite  area,  Ai,  and  strength,  , 
locally  distributed  according  to  a  radially-S3nnmetric  core  function,  fg.  The  discretization 
function  is  also  used  to  reconstruct  the  discretized  quantities  at  later  times.  In  particular. 


^  In  the  case  of  the  shear  layer  flow  the  scaling  values  are  those  of  the  two  free  streams  -  see  Section  IV.  It 
should  also  be  noted  that  when  the  two  scaling  values  of  a  given  property  are  equal  then  its  field 
distribution  is  trivially  reduced  to  a  constant. 
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(12) 


N 

C(x,t)  =  X  Ci(t)  Ai(t)  fgdx-Xvl)  , 
i=l 

where  Xv  =  Xi.(x>t)  is  the  element  location.-  The  core  function,  characterized  by  the  core 
radius  6  within  which  the  most  significant  contribution  of  each  element  is  concentrated,  is 
a  second-order  Gaussian:  fsCr)  =  — i— exp(-  which  leads  to  second-order  accuracy  when 

core  overlap  among  neighboring  elements  is  maintained  [12], 

Discretization  may  follow  one  of  two  approaches.  In  the  first,  g  are  found  from 

Eq.  (12).  This  yields  high  discretization  accuracy  at  the  cost  that  the  element  strengths  do 
not  necessarily  have  a  physical  meaning  but  are  merely  weights  of  an  interpolation 
function.  It  is  not,  therefore,  obvious  that  the  strengths  should  be  governed  by  the  same 
transport  equations  as  the  properties  they  are  used  to  represent.  This  is  particularly  the  case 
when  transport  equations  which  involve  source  terms  are  considered. 

In  the  second  approach,  Cr  =  C(Xi  )•  This  alleviates  the  problem  encountered  in 
the  first  approach  but  at  the  cost  of  reduced  discretization  accuracy^.  It  is  found,  however, 
that  if  the  discretized  field  is  subsequenlty  integrated,  as  is  the  case  for  the  vorticity  and  the 
P -gradient,  the  error  in  the  integrated  quantity,  i.e.  the  velocity  and  the  scalar,  respectively, 
is  much  smaller^ .  For  this  reason,  in  this  work,  this  second  approach  is  followed  for  the 
vorticity  and  the  p -gradient.  For  the  product  mass-fraction,  for  which  no  integration  of  the 
solution  is  needed,  an  approach  similar  to  the  first  is  implemented  but  which  does  not 
involve  the  drawback  mentioned.  Details  of  this  approach  are  given  in  the  section  dealing 
with  the  integration  of  the  transport  equations. 

Calculation  of  the  velocity  and  scalar  fields:  This  essentially  entails  the  solution  of  a 
number  of  Poisson  and  Laplace  equations.  For  the  former.  Green's  function  based 
solutions  are  invoked,  while  for  the  latter  conformal  mapping  techniques  are  implemented. 

The  Green's  function  solution  of  Eq.(5)  for  the  streamfunction  and  the  subcequent 
differenciation  for  the  velocity  are  lumped  into  one  step  to  yield  a  Biot-Savart  convolution. 
Using  the  discrete  vorticity  field  Eq.  (12)  we  get 

N 

Uco(x.t)  =  2  ri(t)  K5(x-x,(t))  ,  (13) 

i=l 


^  The  discretization  accuracy  in  this  case  is  a  strong  function  of  the  overlap  of  the  cores,  generally 
decresasing  as  the  overlap  increases. 

^  A  reduction  of  one  order  of  magnitude  (from  about  5%  to  0.1%)  in  the  maximum  pointwise  discretization 
error  was  experienced  between  a  Gaussian  vorticity  and  the  corresponding  errorfunction  velocity  profiles. 
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where  FiCt)  =  C0i(t)Ai(t)  is  the  element  circulation  KsCx)  =  K(x)  F(^),  where 
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K(x)  =  -^^andF(r) 
27tr2 


=  2jt  I  r  f(r)  dr  = 

Jo 


=  1  -  exp(-  r2) . 


A  similar  approach  is  followed  in  the  solution  of  Eq.(6).  The  calculation  of  the 
expansion  potential,  like  that  of  the  streamfunction,  is  bypassed  via  a  convolution  which 
directly  relates  the  expansion-source  to  the  expansion  velocity.  A  desingularized  equivalent 
of  this  convolution  is  obtained  using  the  assumed  expansion-source  field  Eq.(12 ),  and  is: 


Ue(x,t)  =  -  X  (-  Ai(t)  VG5(X-Xi(t))  . 


(14) 


where 


VGgCx)  =  VG(x)  F(^)  and  VG(x)  =  -  ^ 

8  27tr2 


The  velocity  induced  by  the  boundary  conditions  is  obtained  from  Eq.(7)  using 
Schwartz-Christoffel  conformal  mapping  and  the  method  of  images  (Section  IV). 

The  P  S-Z  scalar  solution  is  obtained  from  the  corresponding  gradient  field  by 
recognizing  the  scalar  and  gradient  can  also  be  related  via  a  Poisson  equation,  V  P  =  V  -  g  . 

The  Green's  function  solution  of  this  equation  is  further  manipulated  by  performing  a  by¬ 
parts  integration,  using  the  divergence  theorem  and  recognizing  that  the  gradient  field 

decays  at  infinity  (unbounded  domain):  p(x)  =  |'  VG(x-x')  g(x')  dx'.  Using  the 

gradient  field  of  Eq.(12)  a  discrete  desingularized  form  is  obtained: 

N 

p(x,t)  =  X  gi(0  Ai(t)  ■  VG5(x-x,(t))  +  Pb  .  (15) 

i=l 

pb  here  is  the  integration  related  function.  It  is  determined  using  a  Shwartz-Christoffel 
conformal  mapping  while  satisfying  the  boundary  condition  using  images. 

Integration  of  the  transport  equations:  The  time  evolution  of  the  flow  and  scalar  fields  is 
established  by  numerically  integrating  the  transport  equations  locally  for  each  element.  For 
this  purpose,  the  integration  takes  place  in  two  fractional  steps  [14],  the  first  deals  with  all 
processes  other  than  diffusion,  while  the  second  concentrates  on  the  latter. 

During  the  first  step  the  element  location  is  updated  by  numerically  integrating 

^  =  u{Xi,t).  (16) 
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using  a  predictor-corrector  scheme.  This  defines  element  trajectories  and  material  lines. 


The  element  vorticity  is  updated  by  integrating  the  circulation  equation: 
dTjCt)  _  [Vp]i(t)  „du(Xi,t)  ^ 


(17) 


dt  p(Xi,t)  dt 

were  the  pressure  gradient  has  been  replaced  by  the  material  acceleration  using  the 
momentum  equation.  The  integration  is  accomplished  using  a  predictor-corrector  scheme 
in  which  the  material  acceleration  is  established  via  a  two-step  iteration,  forward-difference 
algorithm.  Equation  (17)  is  obtained  using  kinematics'^ . 

Using  kinematics  to  simplify  the  transport  is  also  followed  in  the  integration  of  the 
p-gradient  transport  equation.  The  first  fi'actional  step  in  the  integration  of  this  equation 
describes  pure  convection.  Recognizing  that  in  such  a  case  isoscalar  and  material  lines 
coincide  and  considering  the  kinematical  evolution  of  the  latter  in  conjunction  with 
continuity  we  can  rewrite  the  governing  transport  equation*  (left  hand  side  of  Eq.(15))  as 

gi(t) 


dt 


p()Ci.t)  dli(t) 


=  0, 


(18) 


with  gi  =  giiii  and  dli  =  [dlj ,  where  lii  is  the  unit  vector  normal  to  the  material  line  and  dli 

gi(t) 


is  a  material  line  elemental  segment.  Equation  (18)  implies  that 


p(Xi.t)  dli(t) 


=  constant 


along  a  material-isoscalar  line.  The  constant  is  specified  by  the  initial  conditions,  and  dli 
as  well  as  n  i  are  readily  available  due  to  the  Lagrangian  nature  of  the  scheme. 


When  the  reaction  rate  occurs  at  infinite  rate,  integration  of  Eqs.(16-18)  completes 
the  first  fractional  integration  step.  For  the  general,  finite  reaction  rate  case,  however,  the 
integration  of  the  product  mass-fraction  transport  equation  is  also  carried  out.  In  this  case, 
the  concept  of  product  particles  is  introduced.  These  particles  are  located  at  the  element 


^  Material  line  kinematics  expressed  in  the  form  =  dl  •  Vu,  dl  being  a  material  line  elemental 

dt 

segment,  combined  with  continuity  equation  yields  V  u  =  -L  where  dA  =  dl  dn  is  an  elemental  fluid 

dA  dt 

element  area  with  dl  and  dn  being  the  element’s  dimensions  along  and  normal  to  the  material  line, 
respectively.  Substitutiing  this  expression  into  the  vorticity  equation  and  assuming  that  the  vorticity  is 

constant  over  the  area  of  the  fluid  element  yields:  da  +  (V-u)  co  =  -1-  (dA  da  +  d(d^).  (^  |  =  d(.adA)  _  dr 

dt  dA\  dt  dt  /  dt  dt 

*  For  derivation  see  Ref.  [4] 
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centers  and  directly  experience  the  reaction  process.  The  need  to  use  these  particles  stems 
from  the  fact  that  the  nonlinearity  of  the  reaction  requires  the  use  of  precise  values  for  the 
properties  involved  in  the  reaction  rate.  As  already  seen,  the  strengths  of  the  elements  are 
crude  approximations  of  the  values  of  the  field  properties  at  the  same  location. 
Implementing  reaction  on  the  element  strength  will  yield  erroneous,  even  non-physical 
results.  In  addition,  depending  on  the  diffusion  algorithm  to  be  used,  the  strengths  may 
not  experience  the  diffusion  process  (e.g.  in  the  traditional  expanding  core  [13],  the 
random  walk  [1]  or  the  diffusion  velocity  [15]  schemes).  Since  only  material  which  mixes 
molecularly  can  bum,  the  use  of  the  strengths  for  reaction  will  unavoidably  lead  to  errors^ . 

This  last  point  is  important  in  the  implementation  of  diffusion.  Problems  related  to 
the  traditional  core  expansion  scheme  are  associated  with  the  transport  of  expanded  cores 
during  convective  [16].  The  implementation  of  diffusion,  on  the  other  hand,  is  exact. 
With  product  particles,  the  expanded  cores  are  not  convected.  As  will  be  seen,  the  cores 
expand  during  the  diffusion  step,  the  new  particle  values  are  established,  the  particles  arc 
convected,  and  then  the  new  interpolation  function  is  created  using  the  original  size  cores. 
The  net  result  of  this  is  that  the  long  term  effect  of  diffusion  is  to  change  the  strengths  of 
the  elements  rather  than  the  cores.  In  this  respect  the  version  of  the  core  expanding  scheme 
presented  here  bears  similarities  to  the  particle  exchange  scheme  of  diffusion  [17]^°. 

The  integration  of  the  product  mass-fraction  equation,  like  that  of  the  rest  of  the 
transport  equations,  is  carried  out  in  two  fractional  steps.  In  the  first  step  we  integrate: 

w(x,,t),  (19) 

for  each  product  particle  via  an  Euler  predictor-corrector  scheme.  This  yields  the  product 
values  prior  to  diffusion,  denoted  as  Yp(%i,t).  These  values  are  used  to  rediscretize  the 
product  field  and  obtain  the  element  product  strengths  as  noted  earlier,  i.e. 

N 

X  Ypi(t)  Ai(t)  f5(lx-Xil)  =  Yj(xi,t)  (20) 

i=l 

The  solution  of  this  system  of  equations  is  simplified  by  taking  into  account  the  nature  of 
the  core-function  which  implies  that  only  the  closest  neighbors  contribute  to  the  strength  of 
an  element.  Thus,  in  establishing  the  strength  of  a  given  element,  only  elements  which  arc 
r/5  <  fn  (i.e.  exp(-(r/5)  )  >  10’^)  away  from  the  element  are  considered.  A  function 


^  In  this  case,  reaction  will  cease  as  soon  as  any  originally  premixed  material  (i.e.  from  initial  or  inlet 
conditions),  bums. 

10  We  are  indebted  to  Dr.  Georges-Henri  Cottet  for  pointing  this  out  during  a  recent  discussion. 
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iteration  is  used  in  the  solution  of  Eq.(20)  which  is  greatly  helped  by  the  availability  of 
good  initial  guesses,  offered  by  the  values  of  the  strengths  at  the  previous  time  step. 

To  check  the  discretization  error  which  by  construction  is  small,  and  to  minimize  its 
effects,  the  product  particle  values  are  reconstructed  from  the  obtained  strengths  yielding 
Yp*(%i,t)  with  |Yp*(Xi,t)  -  Yp(Xi,t)|  being  the  value  of  pointwise  discretization  error. 

Subsequently,  the  diffusion  step  is  performed  via  core  expansion  scheme 

5?(t+At)  =  6?(t)+4^,  (21) 

where  C  =  Re  for  ^  =  co  and  C  =  Pe  for  ^  =  g  or  Yp,  which  is  the  analytical  solution  of  the 
diffusion  equation  for  each  element.  For  the  product  mass-fraction,  one  more  manipulation 
is  executed  to  reduce  the  impact  (if  any)  of  the  discretization  error. 

Clearly,  the  solution  after  the  diffusion  step,  Yp^'^(Xi,t),  includes  this  error.  By 
subtracting  the  discretization  error  from  this  solution,  a  new,  corrected,  solution  is  obtained 

Yp(Xi,t)  =  Y;(Xi,t)  +  AY5"®(Xi,t)  (22) 

where  AYp^'^(Xi,t)  =  Y^°(%i,t)  -  Yp^(Xi,t)  which  includes  discretization  error  only  in  the 

##D 

change  of  the  product  particle  field  due  to  diffusion  AYp  (Xi,t).  Since  the  physical 
problem  definition  specifies  the  diffusion  to  be  small,  this  approach  allows  even  sigitificant 
discretization  errors  to  be  tolerated.  In  this  work,  however,  this  effect  was  not  exploited 
and  the  average  discretization  error  was  kept  smaller  than  1%.  The  reason  is  while  the 
discretization  error  is  subtracted  from  the  product  particle  solution  (Yp(Xi,t))  it  is  not 
removed  from  the  product  field  solution  (Yp(x,t))  since  it  is  still  included  in  the  product 
strengths.  Thus,  while  the  error  is  not  fed  back  into  the  Lagrangian  calculations  (i.e.  the 
integration  of  the  product  equation)  it  is  present  in  any  instantaneous  product  field. 

The  calculation  of  the  density  gradient  necessary  in  the  integration  of  the  circulation 
equation  implies  knowledge  of  the  product  mass-fraction  gradient  field.  In  the  infinite 
reaction  rate  case  this  is  trivially  established  fi'om  the  SZ  gradients.  For  the  finite  reaction 
rate,  it  is  established  by  differentiating  the  product  mass-fraction  field  given  by  Eq.(18)  as 
suggested  by  Anderson  [18],  i.e. 

N 

VYp(x,t)  =  X  Ypi(t)  Ai(t)  Vf5(x-Xi)  .  (23) 

i=l 

The  severe  stretching  of  the  Lagrangian  mesh  used  in  the  discretization  of  the 
transported  quantities,  which  increases  the  distance  between  neighboring  elements,  may 
lead  to  the  deterioration  of  the  solution  accuracy.  To  overcome  this  problem,  a  scheme  of 
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local  mesh  refinement  is  adopted,  based  on  local  conservation  principles,  whereby  elements 
are  continuously  introduced  and  deleted  to  ensure  core  overlap  [18, 19]. 


IV.  GEOMETRY  AND  BOUNDARY  CONDITIONS 

The  computational  domain  and  some  boundary  conditions  are  shown  in  Fig.l.  The 
shear  layer  evolves  in  a  two-dimensional  channel  of  height  H  and  length  Xmax  >  between 
two  parallel  streams  (1  top,  2  bottom)  which  mix  downstream  of  a  thin  splitter  plate. 
Initially,  both  streams  are  at  the  same  temperature  and  density.  To  po,  and  each  carries  a 
single  reactant.  The  top  and  bottom  walls  are  modeled  as  rigid,  free  slip,  impermeable  and 
adiabatic  planes.  At  the  downstream  section,  a  condition  of  vanishing  vorticity  and  scalar 
gradient  is  used  as  outflow  boundary  condition,  and  is  applied  by  removing  the  elements 
which  cross  the  exit  boundary.  The  error  introduced  by  this  rather  arbitrary  exit  boundary 
condition  was  investigated  by  performing  simulations  with  increasingly  longer  domains 
and  quantitatively  comparing  the  solution  in  their  common  sections.  An  example  from  this 
study  is  shown  in  Fig.2  which  presents  results  of  simulations  of  the  forced  non-exothermic 
flow  for  Xmax  =  7,  top,  and  Xmax  =  5  bottom.  The  shear  layer  is  represented  by  plotting 
the  transport  elements  and  their  local  velocity  vectors  relative  to  the  mean  flow^*  .  We 
found  that  the  difference  between  the  two  solutions  was  negligible  at  about  one  channel 
width  upstream  firom  the  exit  of  the  shorter  domain  (left  dashed  line  in  Fig.2). 


The  inlet  profiles  required  in  the  simulation,  the  vorticity  and  p-gradient,  and,  in  the 


case  of  finite  reaction  rate,  the  product  mass-fraction,  are:  a)(x=0,y,t)  =  - 


_ 


U,-Uo 


fiz  a 


Z(y), 


Vj3(x=0,y,t)  =  ^ Z(y),  Yp(x=0,y,t)  =  Z(y) ,  where  Z(y)  =  exp  . 

a 


In  these  above  expressions,  the  channel  height  is  the  non-dimensionalizing  lengthscale, 

while  Ui,  To  and  po  are  the  velocity,  temperature  and  density  scales,  respectively,  a,  the 

standard  deviation  of  the  Gaussian  profiles,  is  defined  by  requiring  that  two  wavelengths  of 

Y  (x=0  V  t)  =  Y  Z(y) 

the  most  unstable  mode  of  the  uniform  density  shear  layer  fit  p  Pmax  ’ 

within  the  channel  height  (a  =  0.04).  Yj^^  is  chosen  as  0.4. 

The  boundary  conditions  are  satisfied  by  using  a  Shwartz-Christoffel  conformal 
transformation,  which  maps  the  entire  channel  region  onto  the  upper  half  of  a  complex 
plane  (Fig.3).  In  this  mapping  the  two  fluid  streams  appear  as  volume  sources 

The  physical  and  numerical  parameters  used  in  obtaining  the  results  of  Fig.2  not  specified  here  are  the 
same  as  those  of  the  main  body  of  the  results  of  this  paper  which  are  given  in  the  following  sections. 
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S5mQmetrically  located  with  respect  to  the  origin,  and  solutions  to  the  Laplace  equations 
governing  Ub  and  pb  can  be  simply  established  by  using  mass  conservation  arguments. 
The  mapping  however,  is  also  beneficial  to  the  implementation  of  the  image  system  for  the 
transport  elements  which  helps  impose  the  boundary  conditions. 

Initialization  of  the  calculation  is  carried  out  by  assuming  that  the  inlet  conditions 
persist  throughout  the  domain.  Hence,  within  the  support  of  the  vorticity,  p -gradient  and 
product  mass  fraction,  a  square  mesh  of  elements  of  side  h=0.0195  are  distributed  along 
nine  material  layers.  The  value  of  the  core  radius  is  5=0.0234.  It  is  noted  that  the 
discretization  of  the  initial  profiles  may  essentially  be  accomplished  by  the  five  central 
layers.  The  two  extra  layers  on  either  side  of  the  central  five  are  introduced  to 
accommodate  later  rediscredizations  (see  previous  section)  which  may  yield  significant 
values  for  the  strengths  of  the  elements  on  the  extra  layers  as  a  result  of  diffusion.  The 
number  of  extra  layers  to  be  added  is  specified  by  considering  a  diffusion  length  argument 
L  «  VW«  where  tres  =  Xniax/[0.5  (Ui  +  U2)]  is  the  average  residence  time  of  elements 
within  the  computational  domain  and  a  is  the  thermal  (or  any  other)  diffussivity. 

External  forcing  is  implemented  at  the  inlet.  The  forcing  signal  consists  of  in-phase 
components  of  the  most  unstable  mode  of  the  uniform-density  layer  and  its  subharmonic, 
both  at  an  amplitude  Af=As=0.025.  This  forcing  leads  to  early  roll-up  and  pairing  of  the 
uniform  density  layer  [6].  Forcing  is  implemented  by  displacing  elements  at  the  inlet  in  the 
cross-stream  direction  according  to  the  forcing  signal. 


V.  RESULTS  AND  DISCUSSION 

Simulations  using  a  range  of  the  governing  physical  parameters  were  obtained  and 

used  to  analyze  the  flow  and  reacting  field  dynamics  [20-22].  In  another  study  the  effect 

the  inlet  boundary  condition  was  investigated  [23].  In  this  paper  we  concentrate  on  the 

effect  of  the  numerical  parameters  and  we  keep  the  physical  parameters  fixed. 

We  consider  a  shear  layer  with  velocity  ratio  r  =  ^  =  0.5  and  Reynolds  (and 

Ui 

Peeled)  number  Re  =  =  12800,  where  v  is  the  kinematic  viscosity,  evolving  in  a 

channel  of  length  X^ax  =  5.  For  the  finite  reaction  rate  cases  the  normalized  activation 


temperature,  frequency  factor,  enthalpy  of  reaction  and  mass  stoichiometry  ratio  are 
Ta  =  =  10,  Af  =  =  640,  Qo  =  =  6  and  ([)  =  1  respectively.  The  Damkohler 

Riq  U  1  CpTo 

number,  Da  =  -^^,  where  idif  =  —  is  the  diffusion  time  scale  and  Xche  =  ^exp(^)  the 
tche  a  Af  Tf 
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chemical  reaction  time  scale,  and  the  Karlovitz  number,  Ka  =  where  Taw  =  ^  is  the 

Xflw  Ui 


flow  time  scale,  are  Da  »  168100  and  Kad  ®=  0.08,  respectively. 

In  the  simulations  the  spatial  resolution  specified  via  h,  6,  the  number  of  layers 
and  the  parameters  of  the  injection-combination  scheme  was  kept  invariant.  The  temporal 
resolution  was  increasingly  refined  to  test  the  convergence  of  the  solution  in  time.  This 
was  carried  out  for  the  more  general,  finite  reaction  rate  case.  The  time  steps  considered 
were  At  =  0.2,  0.1,  0.05,  0.025.  Figures  4  and  5  present  results  from  this  study.  In  both 
figures  the  At  =  0.025  solution  is  used  as  the  base-comparison  case. 

Figure  4  presents  the  average  pointwise  error  (w.r.t  the  At  =0.025  solution)  of  the 
velocity,  (3  S-Z  variable  and  product  mass  fraction  fields,  i.e.  the  three  variables  related  to 
the  integration  of  the  transport  equations,  for  the  three  different  time-step  cases  considered, 
i.e.  At  =  0.2,  0.1  and  0.05.  In  all  cases  only  the  0  <  x  <  4  section  of  the  computational 
domain  was  considered  as  explained  earlier.  The  solution  used  in  the  calculation  of  the 
error  is  obtained  on  a  9x51  mesh,  as  the  average  of  fifty  instantaneous  solutions  obtained  at 
the  same  times  for  all  numerical  solutions.  The  figure  clearly  shows  that  the  solution 
converges  in  time  and  that  this  convergence  is  approximately  linear. 

More  detail  is  offered  by  Fig.5,  where  the  profiles  for  the  velocity  (top),  the  p  S-Z 
variable  (middle)  and  the  product  mass  fraction  (bottom)  of  the  solution  obtained  as 
described  in  the  previous  paragraph,  are  presented  at  the  same  streamwise  location  x=2. 
Pictures  present  the  comparison  of  the  At  =  0.2  solutions,  with  the  base.  At  =  0.025, 
solution.  The  figure  confirms  the  result  of  the  previous  figure.  We  found  that  the  At  =  0.1 
timestep  yields  solutions  which  essentially  capture  most  of  the  features  of  the  more 
temporally  refined  solution.  As  a  result,  in  what  follows.  At  =  0.1  simulations  will  be  used 
in  the  further  analysis  of  the  flow. 

Figure  6,  displays  two  dimensional  visualizations  of  the  instantaneous  vorticity 
((a)-top)  and  temperature*^  ((b)-top)  for  the  finite  reaction  rate  flow.  These  are  contrasted 
to  corresponding  visualizations  of  the  infinite  reaction  rate  flow  (bottom  of  (a)  and  (b)). 

Substantial  similarity  is  observed  between  the  finite  and  infinite  reaction  rate  results. 
The  finite  reaction  rate  results  suggest  that  the  parameters  chosen  for  the  simulations  are 
such  that  reacting  field  corresponds  to  fast  (as  compared  to  the  flow)  reaction,  i.e.  low 
Karlovitz  number  combustion.  The  infinite  reaction  rate  model  describes  the  extreme  of 
this  case,  i.e.  zero  Karlovitz  number  combustion  in  which  the  flame  is  infinitely  thin  and 


Under  the  assumptions  of  the  physical  model  used  here  the  temperature  fields  are  equivalent  to  those  of 
the  product  mass  fraction. 
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cannot  be  quenched.  Thus,  under  the  conditions  assumed  in  the  simulations  the  infinite 
reaction  rate  model  offers  a  good  approximation  of  the  finite  reaction  rate  flow. 

In  the  context  of  this  discussion  the  similarity  of  the  finite  and  infinite  reaction  rate 
results  is  rriost  important  in  that  it  offers  further  validation  of  the  numerical  methodology. 
From  a  numerical  standpoint  there  are  fundamental  differences  between  the  finite  and 
infinite  reaction  rate  numerical  models  (Section  III).  Most  importantly  are  all  parts  related 
to  the  integration  of  the  product  mass  fraction  equation.  This  equation  is  absent  in  the 
infinite  reaction  rate  model.  Thus,  by  the  simple  fact  of  the  similarity  of  the  finite  and 
infinite  reaction  results  a  qualitative  validation  is  offered  for  the  numerical  integration  of  the 
product  mass  fraction  equation,  the  success  of  the  implementation  of  the  concept  of  product 
particles,  the  new  version  of  the  core  expansion  scheme  which  includes  the  inversion  of 
discretization-interpolation  function,  etc. 


A  more  severe  test  is  demonstrated  in  figure  7  where  the  instantaneous  contribution 


•  H 

of  the  baroclinic  torque,C0B  ~  ^ 


is  plotted  for  both  the  finite  (top)  and  infinite 


(bottom)  reaction  rate  cases.  The  0  <  x  <  3  part  of  the  domain  is  shown.  As  can  be  seen 
from  the  figure,  the  characteristics  of  the  field  contribution  of  the  baroclinic  term  of  the 
vorticity  equation  (Eq.(12))  are  quite  complicated.  This  should  be  expected  as  this  quantity 
represents  the  rate  of  change  of  a  gradient  field.  Satisfactory  explanations  for  it,  however, 
were  proposed  [20].  Our  interest  here  is  the  similarity  between  the  finite  and  infinite 
reaction  rate  results.  The  level  of  detail  in  the  similarity  of  such  a  complicated  field 
substantially  increases  the  confidence  in  the  obtained  solution. 

More  specifically,  the  density  gradient,  required  in  the  calculation  of  cob,  is  related 
to  the  product  mass  fraction  gradient  field  via  the  equation  of  state  and  the  y  S-Z  variable. 
As  was  described  in  Section  III  the  product  mass  fi'action  gradient  field  is  obtained  using 
different  numerical  formulas.  The  fact  that  the  final  result  for  (Ob  is  similar  even  to  the 
smallest  scale,  offers  a  good  consistency  check  between  the  two  approaches. 


VI.  CONCLUSIONS 

The  two-dimensional  Lagrangian  Transport  Element  Method,  an  extension  of  the 
Vortex  Element  Method  able  to  resolve  reacting  scalar  transport,  has  been  presented  and 
implemented  in  the  simulation  of  the  non-premixed  shear  layer  flow.  The  method  is  able  to 
resolve  both  the  variable  density  flow  and  the  associated  exothermic  reacting  field. 

Results  indicate  that  the  method  is  able  to  accurately  reproduce  the  essential  features 
of  the  flow.  Convergence  of  the  solution  in  time  is  presented  and  is  found  to  be 
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approximately  linear.  It  is  shown  that  the  finite  reaction  rate  solutions  obtained  correspond 
to  those  at  low  Karlovitz  number  and  bear  significant  similarity  to  those  of  the  infinite 
reaction  rate  model.  This  similarity  offers  validation  for  the  part  of  the  numerical 
methodology  related  to  the  integration  of  the  product  mass-fraction  equation  since  this  part 
is  not  common  to  the  finite  and  infinite  reaction  rate  models  of  the  reacting  field.  Examples 
of  this  are  the  utilization  of  product  particles  versus  product  elements  for  the  reaction,  the 
use  of  alternative  formulas  for  calculating  the  product  mass-fraction  gradient  and  the 
introduction  of  a  novel  version  of  the  core-expansion  scheme  which  does  not  suffer  from 
the  problems  associated  with  earlier  versions. 
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Xmax 


Figure  1  The  geometry  of  the  computational  domain  together  with  the  initial  element 
configuration  and  some  of  the  boundary  conditions . 
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the  flow  is  forced,  non-exothennic  and  at  t  =  10.  The  length  of  the  domain  is  Xmt 
for  the  top,  and  Xmax  =  5.  for  the  bottom.  The  dashed  lines  are  at  x  =  4. 


Figure  3.  Illustration  of  the  basic  features  of  the  Schwarz-Christoffel  conformal  mapping. 


Solution  Convergence  in  Time 


At 


Figure  4.  Temporal  convergence  of  solution  for  the  forced,  finite  reaction  r^  case  with 
Xniax  =  5.  Convergence  checked  by  comparing  solution  with  the  base 

Ats  =  0.025  case.  The  streamwise  velocity,  P  (or  A,)  S-Z  variable  and  product 

mass  fraction  field  distributions  are  investigated,  i.e.  ^  =  u,  P,  Yp.  Each 
solution  is  established  as  the  average  of  50  instantaneous  solutions  obtained  at 
the  same  times  for  all  cases,  oyer  a  9x51  rectangular  mesh,  i.e.  NP  =  459, 
covering  the0<x<4,  0<y<l  section  of  the  computational  domain. 
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The  product  mass-fraction  (or  temperature)  fields  corresponding  to  the  vorticity  fields  above. 

Figure  6 


Figure  7.  The  baroclinic  vorticity  generation  rate  of  the  reacting  at  finite  (top)  and  inflnite  (bottoin) 
reaction  rate  shear  layers  at  time  t  =  10.  The  0  <  x  <  3  part  of  the  computational  domain 
is  shown.  Solid/dashed  contours  indicate  negative/positive  vorticity  at  an  increment  of  0.4. 


